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STANISŁAW BAŃKA, CEZARY MOSKWA 


INVESTIGATION OF MULTIPURPOSE 
CONTROL SYSTEMS USING THE ASMCS PACKAGE 
FOR ANALYSIS AND SYNTHESIS 
OF MULTIVARIABLE CONTROL SYSTEMS* 


Institute of Control Engineering, Technical University of Szczecin, 
ul. Gen. Sikorskiego 37, 70-313 Szczecin 


Summary. This paper concerns the properties of composite (multipurpose, 
multiloop) control systems for linear multi-input multi-output (MIMO) plants de- 
scribed by square proper rational transfer matrices. In the stochastic case multipur- 
pose control systems simultaneously ensure: (i) complete dynamic decoupling of 
control loops, (ii) arbitrary closed-loop pole placement and (iii) steady state output 
rejection of non-zero mean values of disturbances, as well as zero steady-state regula- 
tion, or tracking errors. Three kinds of control system structure are considered, 
assuming that the state vector x of the plant is either accessible or unaccessible for 
direct measurement. In the latter case, a (stationary) Kalman filter is applied to 
obtain an optimal estimate of the plant's state and output vectors £ and j, which are 
then used in state variable feedback and to define the error signal vectors e = yo—5 
and e = yę—y respectively. Numerical example and results of simulations are in- 
cluded to illustrate how multipurpose control systems work and to examine proper- 
ties of the proposed control system structures from a practical points of view. 


Key words: multivariable control systems — dynamical decoupling — pole 
placement — state estimation — polynomial matrix approach. 


INTRODUCTION 


The first papers on multipurpose control systems appeared in the 80's. 
The majority of them concerned optimal (LQ, LQG) control systems in which 
weighting matrices of their LSQ control cost functions were suitably chosen to 
produce desired pole location of the closed-loop control systems (see for ex- 
ample [1, 11, 18]). 

The next steps in this direction were made in the papers [10] and [20], 
where the above-mentioned weighting matrices were chosen to achieve, 


* This work has been partially supported by the Committee for Scientific Research of the 
Republic of Poland under Grant No 8850502205. 


simultaneously, input-output decoupling, complete and arbitrary closed-loop 
pole placement with disturbance rejection and reference signal tracking. Unfor- 
tunately, it has been proved by Hautus and Van der Wegen [13], that there are 
cases when it is not possible to fulfil all of these requirements with such 
a simple (single-loop) control system structure. 

The first results, in which the above purposes were truly achieved, were 
obtained by Wolovich [21]. This paper was concerned with square, non-sin- 
gular, multi-input, multi-output, strictly proper plants with prescribed class 
of non-decreasing deterministic inputs (including disturbance and reference 
signals). 

Many attempts have been made by Bańka, for example [2, 3, 6], to extend 
the results of Wolovich to more general, proper plants with stochastic and 
deterministic disturbances. Bańka, [3, 6], has shown that all of the design goals 
can be achieved also for nonsquare proper plants in the stochastic and deter- 
ministic case. However, properties of the designed control systems from prac- 
tical point of view with reference to the more complex (multi-loop, multi-part) 
control system structure are poorly understood even those known for square 
plants. 

In this paper we are going to analyze and illustrate how the proposed 
multipurpose control systems work in real situations for the class of proper 
(square) continuous plants with stochastic and deterministic disturbances, with 
both regulation and tracking tasks. 

The paper considers three cases of a multipurpose control system structure 
as follows: 

(1) multivariable control system structure with access to the original output 
and accessible state vector of the plant; 

(2) multivariable control system with filtered output and filtered state vec- 
tor obtained from a stationary Kalman filter; 

(3) multivariable control system structure using the original output of the 
plant and the estimated state vector. 

In all cases it is assumed that the multipurpose control systems should 
simultaneously achieve: 

(i) complete dynamic decoupling of the closed-loop control system, 

(ii) arbitrary closed-loop pole placement with internal stability and internal 
property of the control system, and 

(iii) complete steady state output rejection of non-zero mean values 
of stochastic disturbances and zero steady state regulation or tracking 
errors. 

In order to solve this problem we apply the direct polynomial matrix ap- 
proach using the algorithm given by Bańka in [3] (see also [2, 21]). In this 
algorithm the following principles are applied simultaneously: diagonal dyna- 
mic decoupling and pole placement by state variable feedback with „input 
dynamics” [22], the „internal model” principle in the form of sufficient 
conditions given by Callier and Desoer [7] and the well-known „separation” 
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principle to solve the optimal (steady-state) estimation problem irrespective of 
solving the control tasks. 

Investigations are illustrated by design and simulations for three cases of 
control systems for two-input, two-output second order linear plant using the 
ASMCS package v. 2.03 — 31.12.1992 developed by the authors for analysis and 
synthesis of multivariable control systems. 


THE PLANT AND STRUCTURE OF THE CONTROL SYSTEMS 


We consider a controllable and observable linear multivariable plant 
(which may be unstable, non-minimum phase or both) defined by the state and 
output equations: 


X(t) = Ax(t) + Bu(t)+ Er(t) + Gwt(t), (1) 
y(t) = Cx(t)+ Du(t) + v(t), (2) 
where x(t)e%"”, u(t)je%" and y(t)e%'(m =1) are the state, input and output 


vectors respectively. The vectors w(t)e R?” and v(t)e%' represent zero-mean 
white random Gaussian processes with the mutual covariance matrix 


z|| „4 be" et] 7 Ę y |50-0 (3) 


which can be semipositive definite for the positive-definite density matrix 
V > 0. Non-zero (non-diminishing) mean values of the disturbances w(t) are 
represented by the vector r(t)e 7”. Whenever E = G the latter can also re- 
present additional deterministic disturbance signals r(t)eR' (r ź£ p). 

Since the polynomial matrix approach in the frequency s-domain is adop- 
ted, the plant model is transformed to 


y = B,(s)A; '(s)u+ Az '(5)B,(s)w +43 *'(5)B, (S)F+v, (4) 
where the proper rational transfer matrix 
B,(5)Ar '(5) = C(sl,—4)71B+D (5) 


between the signals u and y is defined by relatively right prime (RRP) 
polynomial matrices B,(s)e%Ż[s]'"" and A,(s)e%[s]”*”, with deg,,[B,(s)] 
< deg.,[A,(s)], i= 1, 2,..., m and 4A,(s) column-reduced. 

AII remaining transfer matrices in the eq. (4) are strictly proper and defined 
by the relatively left prime (RLP) matrix fractions descriptions: 


A> '(s)B,(s) = C(sl,—A) 'G (6) 
and 
Az '(5)B>(s) = C(sI,—4) 'Er(s), (7) 


where A,(s)ef8[s]'*' and B,(s)e%[s]'*” are RLP with A,(s) row-reduced as 
well as A;(sje2[s]'"' and B,(sje%[s]'*"” are RLP with 4A;(s) row-reduced. 
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The column structure of the denominator matrix A,(s) is determined by its 
column degrees deg,,[4,(s)] = £,, i = 1, 2,..., m, and by the (non-singular) 
highest-column-degree coefficient matrix I'_(A, (s)). Both.of these depend on the 
controllable pair (A, B) of the state equation (1). By analogy, the row structure 
of the matrix 4,(s) is defined by its row degrees deg,,[4,(5S)] = v,, 
i= |, 2,..., l, and the highest-row-degree coefficient matrix I„(4, (s)). These 
depend on the observable pair (A, C) of the plant description. All the assumed 
properties of the above-mentioned polynomial and transfer matrices are satis- 
fied according to the well-known Wolovich structure theorem [16, 22]. 

Note that the transformed disturbance vector r(s) is included in the transfer 
function matrix (7), and so the vector r in (4) is a „fictitious” (impulsive) input 
signal for the deterministic disturbance model. For tracking problems it is 
assumed that the reference signal vector y,(s) is generated from a reference 
model defined by 

Vos) = 4o '()I,jo, (8) 


where Ag '(s)I,, is a diagonal and strictly proper transfer function matrix with 
the impulsive input signal y,. We define polynomials m, (s), mą4(s) and 
my + (s) as the completely unstable least common denominators of the unstable 
parts of the transfer matrices (6), (7) and (8) respectively. 

The three possible structures of the control system described in the time 
domain, are presented in Figs. 1, 2 and 3. In Fig. 1, where direct state feedback 
f =Fx is used, the proposed multipurpose „two-part” compensator/controller 
includes a diagonal controller and a feedforward compensator („input dyna- 
mics”). If the plant's state vector is not available for direct measurement then 
the control system structures, pictured in Figs. 2 and 3 include three parts 
which are: a diagonal controller, a feedforward compensator and a stationary 
Kalman filter. The Kalman filter in Fig. 2 has two-input signal vectors u and y, 
and two-output signal vectors f and $. In Fig. 3 the filter has only one output, 


which is the feedback signal vector f =FX. 
DYNAMIC 


Fig. 1. Structure of the control system with accessible state vector of the plant 


11 


THE PLANT 


Fig. 2. Structure of the control system with Kalman filter and feedback from filtered output of the 
plant 


DYNAMIC 
FEEDFORWARD 


CONTROLLER 


KALMAN 
FILTER 


Fig. 3. Structure of the control system with Kalman filter and direct feedback from output of the 
plant 


The signal vector f = FX (and f =Fx in Fig. 1, respectively) is employed 
primarily (but not only) to decouple the „inner” part of the control system 
between the signals q and y, and in the case of Fig. 2 between q and y. The 
estimated plant output vector ) = C£ + Du, in Fig. 2, is used to define the error 
signal vector e = yg—y. In Fig. 3 the error signal vector is defined simply as 
AGA ANIE A, 

In cases when the plant's state vector x is unaccessible, the estimated state 


A 


vector £ is obtained by minimizing the LSQ cost functional 
Ą l j 00 
J = wE[[x0—+(01C:0)-+(017) = > |-jo tr7/wu(s)ds. (9) 
where ) ,,(s) is the spectral density matrix of the state estimation errors 
e,=x—-X. 
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The most general structure of the multipurpose compensator/controller, 
given by the time domain descriptions depicted in Figs. 1 — 3, is furnished in the 
frequency s-domain in Fig. 4. It is suitably defined by: 

— the strictly proper transfer matrix M; '(s)N>(s) for the diagonal control- 
ler, 

— the proper and low-order transfer matrix G' '(s)L(s) for the dynamic 
feedforward compensator („input dynamics ”), and for the case when the plant's 
state vector is unaccessible by: 

— the strictly proper transfer matrices 


Q"'(s)H(s) = F(sl,—A+KC) 'K (10) 
and 
O '(s)K(s) = K(sl,—A+KC) '(B-KD) (11) 
for the Kalman filter with the feedback matrix F, 


THE PLANT | 
I B (s) | 
1 MODEL OF REFERENCE 3 
| SIGNALS * +: | 
B,(s) [— | s) 
-] ma p IE "3 ż 
5 i 


CONTROLLER DYNAMIC FEEDFORWARD 


| | | |o Ą 


BI 
Q(s) 
KALMAN FILTER 


Fig. 4. The full structure of the control system in s-domain 
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— and, in the case of Fig. 2 only, by the strictly proper transfer matrix 
Q; '(s)H,(s) = C(sl,—A+KC) 'K (12) 

and proper transfer matrix 
Q; '(s)K, (s) = C(sl,— A+KC) '(B-KD)+D (13) 


for the Kalman filter along with the plants output matrix C. 

The main problem is stated as follows. Given the transfer matrices of 
a plant, defined by the equations (4)—(7) and (8), determine the transfer mat- 
rices Mz *(s)N2(5), G”'(5)L(s), Q”'(s)[H(s)|K(s)] and Q; '(s)[H,(s)|K,(s)] as 
components of the multipurpose compensator/controller so as to realize the 
desired design objectives (i) — (iii). 

AII the required transfer matrices are calculated in a RLP-factored 
form with non-singular row-reduced denominator matrices M,(sje[s]'*', 
G(s)eR[s]”*"”, O(sjeR[s]"*"”" and Q,(sje%[s]'*"'. The latter denominator 
matrices for the Kalman filter must be (Hurwitz) stable with the same zeros of 
their determinants |Q(s)| = «|Q,(s)|. 

Row structure of the polynomial matrices Q(sjeM[s]" ”" and Q,(s) 
eR[s]'"' usually differs from the row structure of the denominator matrix 
A,(sjeM[s]'"' even for the case of I = m. They depend on the (observable) 
pairs (F, A KC) and (C, A—KC) respectively. Since they cannot generally be 
regular, minimal degree solutions of appropriate polynomial matrix equations 
may not be unique [8, 12, 14—17, 19]. Therefore, special efforts have been 
made to avoid necessity of solving them [3, 5]. The algorithm used in this 
paper quarantees that the transfer matrices (10) —(13) are unique and strictly 
proper without solving any polynomial equation (apart from the transfer mat- 
rix (13) which is always proper for proper plants). 

We assume that the polynomial matrices M,(s) and L(s) as well as M,(s) 
and Q(s) are RRP. 


LAYOUT OF THE DESIGN METHOD 


The method and the design algorithms for solving the multipurpose control 
problem described above, with square plants have been presented in [2, 3, 4]. 
Details of, these are omitted. Loosely speaking, the idea of the method is as 
follows. Using the linear state variable feedback, along with the dynamic feed- 
forward G "(s)L(s), we decouple the „inner” part of the system between the 
signals q and y to obtain the diagonal transfer matrix 


T,„(s) = B,(s)[G(s)A,(s)—F(s)] "' Ł(s) = N(s)D”"'(s) (14) 
where 
D''(s)=J"'(s)[G(s)4,(5)—F(s)] ' L(s) = 
= [(6(8)4,6)-F(s))J(s)] "" L(s) (15) 
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and 


J(S)N”"(s) = By '(8) (16) 


with determinants |J(s)| = |L(s)|. Since zeros of these determinants are equal 
to the transmission zeros of the plant defined by polynomial matrix B,(s), 
we have cancellation of unobservable and uncontrollable zeros in „inner” 
decoupled part of the system. Note, if the plant is non-minimum phase then the 
control system does not satisfy the requirement of internal stability. 

Next we apply the „internal model” principle to design a feedback system 
with the „decoupled plant” defined by the eq. (14) and the diagonal controller 
M; '(s)N,(s). According to the sufficient conditions of this principle, given by 
Callier and Desoer [7], the denominator matrix of the controller is chosen as 
M, (s) = diag [m(s)] (see Wolovich [21]), where m(s) is the least common multi- 
plier of the polynomials m++(2), m3,(s) and mo9+(s). Hence, we obtain the 
diagonal Diophantine polynomial matrix equation 


M,(s)D(s)+ N(s)N; (s) = 4(5), (17) 


for given M,(s), N(s) and 4(s), where elements of the 4(s) = diag[0,(s)], 
i=1,2,...,ł, are chosen to be stable (monic) polynomials of degree 
deg [m(s)] + deg [d,(s)] matched to the assumed configurations of the clo- 
sed-loop control system poles. The minimal degree solution of the eq. (17) 
simultaneously yields both the numerator matrix N,(s) for the controller 
and the denominator matrix D(s) for the decoupled „inner” part of the 
system. 

The design algorithm developed by using this idea is presented in the works 
[2] and [3] (see also [21 ]). This algorithm forms part of a software package for 
synthesis and simulations of multivariable control systems (the last release is 
ASMCS v. 2.03 — 31.12.1992) developed by the authors, which is implementable 
on IBM PC with Hercules and VGA/EGA cards. The package contains many 
subroutines for calculating MFDs of linear MIMO systems and for solving 
Diophantine matrix polynomial equations and spectral factorization problems, 
as well as the subroutines needed for operations on polynomial matrices. Also, 
it contains the profesional software for simulating and editing the data and 
graphic results obtained. The maximum size of jobs realized by ASMCS is 
restricted to linear plants of up to nine degrees with six inputs and six 
outputs. 

Dynamical behaviour of the overall multipurpose control system (in the 
general case when the plant state vector is unaccessible) is described in the 
s-domain by the block polynomial matrix forms 


P(s)X, = Q(U, (18) 


Y =R(s)X,+V(s)U, (19) 
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where * 
GA, —L —l[, 0 0 0 0 
0 M, 0 —N, 0 0 0 
—KA, -HB;, 0 (0) 0 —H -H 0 
P(s) =| -K,A,-H,B, 0 0 —Q, -H, -H, 0 
0 0 0 0 A, 0 0 
0 0 0 0 0 A 0 
0 0 0 0 0 0 A; 
0.0 0 0 x, 
N, 0 0 0 'q 3 
0.00 0 H f * e 
Q(s)=|0 0 0 H,|], X=|5 | U= "ak Y=|y|, 
0 0 B, 0 w i 5 
0 B, 0 0 r 
I, 0 0 0 Jo 
and 
0.00 77, 0 0 1 000 0 
Rs)=|B, 0 0 0 I, I, 0 and V(s)=|0 0 0 4,|. 
0.00 7, 0 0 0 000 0 


Treating this description as a left-right polynomial matrix fraction of the 
generalized (block) rational transfer matrix between the block vectors Y and 
U defined by 


<= 
= 


0 


e; 
T(s) =R()P"'()0(5)+V(s) =| 7; (20) 
K; 


AAA 
RAR 
BA 


=! 
= 


we obtain the following transfer matrices: 
Ty, = I,—Q; '(K,A,+H,B,)4; 'QLN,AG"'L,, (21) 
T Ę” [Q: '(K,A, + H,B,)4; '(QLN,Q; 'H,—HM,) 
—Q,'H,]4A;'B;, (22) 
T„= [Q; '(K,A, +H,B,)4; '(QLN,Q;*H,—HM)) 
—Q,'H,]4;'B,, (23) 


* [he operator s is ommitted for the interests of simplicity. 
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T„ =[Q; '(K,A,+H,B,)4;'(QLN,Q;'H,—HM,)-Q;'H,], (24) 


T,y, = B,Ar"QLN,AG 'L,, (25) 
T,, = [I, —B, 47 '(QLN,Q; 'H, —HM,)]A; 'B;, (26) 
T,, = [I,—B, 47 '(QLN,Qr'H, —HM,)]A; 'B,, (27) 
T,, = [I,—B, 4; '(QLN,Q; 'H,—HM])], (28) 
Ty, = Qr '(K,A,+H,B,)4; 'QLN,A6 '1,, (29) 
and 

EB; 4 Tw; By, a | Re and Ę, " (A (30) 

where 


A,(s) = m(s)(QGA, -KA,-HB,)+QLN,Q;'(K,A,+H,B,). (31) 


In the general case of the considered problem, with unaccessible plant's 
state vector, it has been shown in [2], [3] and [5] that the polynomial matrices 
O(s), H(s) and K(5) of transfer matrices defined by the eqs. (10) and (11) for the 
Kalman filter, should satisfy the matrix polynomial equation 


K(s)A,(s)+H(s)B,(s) = O(s)F(5). (32) 
Similarly, the matrices Q,(s), H,(s) and K,(s) in the eqs. (12) and (13) 
should satisfy another matrix polynomial equation 


K,(s)4,(5)+H,(S)B; (s) = Q,(5)B, (3). (33) 


If a solution of the state estimation problem exists then the zeros of the 
determinants 


I26)| = alQ,(5)| = BlsI,—A+KC| (34) 


lie in the stable half-plane Re(s) < 0. Moreover, we have shown also in [2] and 
[3] that 


A,(5) = ©6)L(S)4A(s)J" (s) (35) 
and since |L(s)| = |J(s) we have 
14,(8)| =|2(9)|-|4(5)|. (36) 


In the last, |4(s)| is an arbitrarily assignable Hurwitz polynomial. However, 
cancellation of zeros of |L(s)| and |J(s)| in the eqs. (35) and (15) can be unstable, 
if the plant is non-minimum phase. This fact can cause instability of the overall 
closed-loop control system. 

Except for this latter case, after suitable (stable) cancellations in the eqs. 
(21) —(30) and (31) we obtain, among others, the following diagonal and stable 
transfer matrices 


[6 = 7,6) = N()47-(9N2(5), (37) 
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of the closed-loop system between the reference signals yy and the output 
signals y and the filtered output signals j. Using the Diophantine equation 
(17) it follows that the poles of these matrices are given by the zeros 
of |4(s)|. 

Zeros of the closed-loop control system are determined by the zeros of 
|N(s)|, which consist of plant's transmission zeros and by the zeros of |N,(S)|, 
which are zeros of the controller. Note, the latter can result in a non-minimum 
phase control system even if the plant is minimum phase. 

Since all „internal” poles of the control system defined by the determi- 
nants 


I21(5-I4r(8)] = IQ1(91-I2(5)|-|4(8)1, 


are stable and all of unstable modes of „external” models defined by the eqs. (6), 
(7) and (8) are cancelled with the zeros of m(S) resulting in stable and strictly 
proper transfer matrices defined by the eqs. (21)—(24), we have 


E(lim e(t)) = 0. 


t>o 


This implies zero steady state regulation or tracking errors as well as complete 
steady state output rejection of non-decreasing disturbance mean values defi- 
ned by the vector r(t). Since the systems are also decoupled, we have therefore 
achieved all the design goals (i) —(iii) for any proper (square) but minimum 
phase plant. 


THE RESULTS OF EXAMINATION 


There are many issues related to the design of multipurpose control sys- 
tems. In this paper we restrict ourselves to these issues which are important 
from a practical point of view as follows: 

Q1. When should the full control system structure depicted in Figs. 2 and 
4 be used and what are resulting benefits? 

Q2. What worsens the quality of control processes, in comparison to the 
control system structures with direct feedback from either the plant's state and 
output vectors or both? What is the cause of this worsening? 

Q3. In which cases can the solutions obtained be considered unsatisfactory 
from a practical point of view? 

Q4. How sensitive is the control system design to inaccurate knowledge of 
plant parameters? 

Q5. Is the control system really unstable if a (square) plant is non-mini- 
mum phase? 

In order to address these issues we will design and simulate the a- 
bove-mentioned control structures for a tracking problem with presence of 
stochastic and deterministic disturbances, for the following example. 


2 — Szczecińskie Roczniki Naukowe t. IX, z. 1 
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The Example 


Let us consider two-input, two-output, proper (stable and minimum phase) 
plant described according to the eqs. (1)—(2) by the matrices 


a-|5os] e-[o] olo] 2-[121] 


and 


with white gaussian noises w(t) and v(t) described by the positive-definite mutu- 
al covariance matrix 


1010 0 
Wasi Ssa SE m 
ST p 0.0101 0 i 

0.010 01 


Additionally, the signals r,(t) = 5sin(2n/24t) and r,(t) = —5sin(2n/24t) are 
assumed as components of the deterministic disturbance vector r(t). The de- 
signed control system should simultaneously track the ramp reference signals 
yo(t) and reject above-mentioned components of the deterministic disturbance 
vector r(t). 

Assuming the same closed-loop pole location (i.e. zeros of |4(s)|) for two 
input-output loops of the system as follows: s, = —1, s, = —2, Sz = —3, 
S4 = —4, 85 = —5and s, = —6, we obtain by use of the design algorithm the 
following results: 

— the diagonal matrix transfer function for the controller 


16.305s* + 87.405s* + 151.56s + 80 
s* —0.068539s* j 


M; '(5)N, (s) = diag | 


— the proper transfer matrix for the dynamic feedforward compensator 


s+ 2.6953 — IB pz | 


—1 A 
4 (9L69 = | 3  <s+46058] |—S :-3 


— and the feedback matrix 


_ [2.2504 —0.3047 
" [0.9141 1.9457 | 


As a result of calculations we have designed completely decoupled stable 
and minimum phase closed-loop system with zeros defined by zeros of |N(s)|, 
which are transmission zeros of the plant equal sf > = —2.5+j1.6583, and 
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the additional zeros s34 = —2.1890+j0.44854 and s$ = —0.9827 caused by 
zeros of |N,(s)| for the controller. In this example we have obtained the same 
transfer functions for each input-output loops defined by the eq. (37) since the 
same configurations of closed-loop poles were assured. IDecoupled „internal” 
part of the system between the signals q(t) and y(t) has also been stable (and 
minimum phase) with poles defined by the values —2.34765 + j0.79296. 
The results of simulation of the designed multipurpose control system are- 
presented in Figs. 5—16. In the case of system structure from Fig. 1, with 
accessible state vector x(t) of the plant, we have obtained a very good tracking 
result which is depicted in Fig. 5. In this figure we show the simulated processes 


12G. 
1GG. | 
8G. 
GQ. 
4G. 


26. i 


—26. A i * : A 
a. 5. 1G. 15. 28. 25. 38. 


Fig. 5. Output signals y and reference signals yo compared in the case of accessible plant 
state vector 


for the output signals y(t) which ideally follow the introduced changes in the 
reference signals y,(t). In Fig. 6 the components of error vector e(t) = yo(t) 
—y(t) against a background of r(t) are presented. 

For the case when the plant's state vector was unaccessible, we additionally 
design the stationary Kalman filter with the gain matrix obtained as 


_ [70.6909 — 1.0622 
-|1.7033 1.0123 |- 


According to the assumed intensities of noises w(t) and v(t) defined by (3), the 
calculated (optimal-in-steady-state) gain matrix K produces the transfer mat- 
rices (10) —(13) of the filter with (uncontrollable) stable poles equal —4.3889 
+j0.8732. Since these are located approximately near to the assumed values 
for the system closed-loop poles, the obtained Kalman filter is not „sufficiently 
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G. 5. 16. 15. 2G. 25. 38. 


Fig. 6. Error signals e = yy —y against a background of disturbance r obtained in the case of 
accessible plant state vector 


fast” to the problem considered. It causes slightly worsening the tracking 
processes, especially when the most general structure of the control system 
from Figs. 2 and 4 is applied. They are shown in the succeeding Figs. 7—11. 


Yo.y 
k 6-17 JB ORO OE EE ZP PCE RENE TONE POREEOBEYCO 


1GG.ł--'--:-«:-:::: R EPEE CZTERY TZN PY M WE WNE CA oazę PETTER AW OEŃ NY 0 ra AŁRZZE AA a ; 


<= OOLANIS M. FE ZW. an > KANE: 
i RODNE i 
GG. ONE: PTOP PPE, ; SR szAA SW" SEAN ę 3 ; 510 4 AR GT 0 . e4E w wa wi sł dą w.ś fofninż 6194-19 44 ay 


Lr! SZAREK RONA BO Lekż PA PAŹ „zB OB 

i (o i ; yozy2 | | 

2. CE SODY AE On NE : 
a PROZAK KRKA, s iWpęraa Pa rasy a ia mi mannaa a awa na wkak kad 

a. 5. 18. 15. 28. 25. 3G. 


Fig. 7. Output signals y and reference signals y, compared in the case of unaccessible plant state 
vector with e = po—y 
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UODO EEEE EEEESZTNKZKENNNEM 


G. 5. 1G. 15. 26. 25. 38. 


Fig. 8. Error signals e = yy—y against a background of disturbance r obtained in the case of 
unaccessible plant state vector 


In Figs. 7 and 8 the output signals y(t) compared with the reference signals 
Yo(t) and the error signals e(t) = yo(t)—y(t) against a background of r(t) are 
presented. They concern the case of Fig. 3 control system structure with f = F£. 
As we can see changes between the graphs in Figs. 5—8 are not remarkable. 


120.807 oda Ks OREW AR PCA EEA ODA AAAA TAA NEBĄ O NoE M Ż 
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7. [0 SPP OP ALE Ri PMA M | EGER Er" 7 BELI NIP EDE Wę ik” 
GE Tatan z: ; Anlóyś Wa -v WATPE IE NSZAEEA 4 Z SORJA EAN EO Ebook, 
2 Volt] 
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Fig. 9. Output signals y and reference signals yy compared in the case of unaccessible plant state 
vector and e = yo—5 


0 
IN0) 


r>(t) 


G. 5. 1G. 15. 28. 25. 38. 


Fig. 10. Error signals e = py —) against a background of disturbance r in the case of unaccessible 
plant state vector 


Similarly, in Figs. 9 and 10 we have shown the tracking processes for the 
same as before reference signals y,(t) in the case when the full control system 
structure from Figs. 2 and 4 is used with f (t) = FX(t) and e(t) = yo(t)— S(t). 
However, in these figures changes are visible. In this case dynamical properties 
of the Kalman filter used in the full control system structure worsened the 
tracking processes. 

In order to explain how dynamical properties of the Kalman filter influence 
the control processes, a comparison of components for actual and estimated 
state vectors x(t) and <X(t) obtained during the simulations are presented in 
Figs. 1l and 12. The last figure concerns the results obtained for grown inten- 
sities of noises w(t) to the values 


For these values we have calculated the new gain matrix 


_| 49166 —8.5037 
_|10.1010 5.1840 |” 
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-4G. 
a. 5. 1G. 15. 20. 25. 38. 


Fig. 11. State vector components of the plant and Kalman filter compared in the case of low 
intensities of noises w 


Fig. 12. State vector components of the plant and Kalman filter compared in the case of grown 
intensities of 'noises w 
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from which yields a „more fast” Kalman filter with poles given by the values 
—8.3685 and — 20.4198. This improves the tracking and regulation processes 
in the case of full control system structure employed. These are depicted in 
Figs. 13 and 14. 


12G. 
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Fig. 13. Output signals y and reference signals yo compared in the case of grown intensity of noises 
w and e=ypqo—j) 


Fig. 14. Error signals e = y—y against a background of disturbance r in the case of grown 
intensities of noises w 
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As far as robustness problems with respect to inaccurate knowledge of 
plant parameters are concerned, they are difficult to fully analyze due to com- 
plexity of the proposed control system structures. Hence, simulations of desig- 
ned systems play the important role. However, for this purpose many parame- 
ters can be diverged in the control systems considered. In this paper we restrict 
ourselves only to investigation of influence of time delays not anticipated in 
control loops of the system. For instance, in the considered example the in- 
fluence of time delay introduced to the control signals we have examined. The 
results of simulation with delayed input signals either u,(t— Tę) or u,(t — T5) 
with 7, = ls, are presented in Figs. 15 and 16 respectively. They have been 
performed with grown intensity of the noises w(t) only. These simulations show 
that for this considerable time delay introduced in the control loops, the system 
remained stable and decoupled with slightly worsening tracking processes for 
the output signals y(t). 

AII of the simulations were started from non-zero initial values of the state 
vector xy = [10, —10]” and zero initial values of state vectors for the Kalman 
filter, controller and feedforward compensator. 

In other examples not presented here, the obtained results are not so good. 
They happen if poles of the Kalman filter are considerably less than the poles 
of closed-loop system. In these examples the structure of control system from 
Figs. 2 and 4 can give completely wrong control processes. Indeed filtered 
outputs y(t) tightly track reference signals y,(t) in opposition to true values of 
the output signals y(t) that are far from them. On the other hand, during 
designing non-minimum phase zeros of the closed-loop system caused 
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Fig. 15. Example of tracking processes obtained for delayed input signal u, (t— Ty) 


Fig. 16. Example of tracking processes obtained for delayed input signal u,(t— Ty) 


by |N„(s)| can occur. Such systems cannot track any reference signals in enough 
satisfied way. To avoid that problem, the zeros of |4(s)| assumed in the eq. (17) 
should be suitably changed. 

On the basis of our experience we have stated the following remarks: 

RI. Note that in all presented figures the tracking processes for particular 
components of y(t) are practically independent of the tracking processes in 
each of other control loops. Therefore, a decoupling property is practically 
irrespective to the used control system structure. It confirms that the proposed 
multipurpose control system designs are generally robust with respect to de- 
coupling and stability requirements. 

R2. The full control system structure depicted in the Figs. 2 and 4 can be 
used only if the obtained Kalman filter is sufficiently fast, i.e. if poles of the 
transfer matrices (10) —(13) obtained are considerably greater than the assumed 
poles of closed-loop control system. It depends on the intensity of noises w(t) 
and v(t). If this condition is not satisfied the structure presented in Fig. 3 (or 
Fig. 1, if it is possible) should be used. The resulting benefits are mainly in 
additional reduction of modulations of the control signals u(t). 

R3. The solutions obtained can be unsatisfactory from practical point of 
view in non-minimum phase zeros of the designed control systems caused by 
|N„(s)| occur. They can be eliminated during solving the Diophantine equation 
(17) by suitable selecting of closed-loop poles of the system. However, some 
restrictions on these selections exist since zeros of the diagonal matrix 4(s) 
should be disjoint with the zeros of the matrix N(s). This condition is sufficient 
but not necessary [4]. 


PAŹ) 


R4. As far as the question of Q5 is concerned, we have proved for many 
examples (omitted here because of lack of the space) that such designed systems 
(with square plants) were never stable if the plant had non-minimum phase 
zeros. "The last problem motivated us to look for a new algorithm for solving 
the considered problem without cancellations of any unstable hidden modes, 
especially in the case of non-minimum phase plants. This has been done by 
Bańka, [3, 6], in the case of proper, non-square (right-invertible) plants. 


CONCLUSIONS 


In the paper we have examined properties of multipurpose control systems 
from a practical point of view, for square linear proper plants. The simulations 
proved that the considered control system structures can achieve all of the 
desired design goals and the designed systems are practically robust, especially 
with respect to the decoupling and stability requirements. Hence, they can be 
employed in real control situations. 

However, the design algorithm used in this paper ensures internal stability 
of these control systems only for non-minimum phase (controllable and observ- 
able) plants, which can be unstable or not. Designs possible to obtain (from 
a theoretical point of view) for non-minimum phase plants are useless because 
they are unstable in practice. 

The new method presented in [6] enables us to make synthesis of multipur- 
pose control systems for the more general case of plants with the number of 
inputs greater than the number of their outputs. This new algorithm ensures 
internal stability and internal property of the control system designed for each 
of the conditions with the both unstable and non-minimum phase plants. 


STANISŁAW BAŃKA, CEZARY MOSKWA 


BADANIA WIELOZADANIOWYCH UKŁADÓW STEROWANIA 
Z WYKORZYSTANIEM PAKIETU ASWUD 
DO ANALIZY I SYNTEZY 
WIELOWYMIAROWYCH UKŁADÓW STEROWANIA 


Streszczenie 


W pracy przedstawiono wyniki badań złożonych, wielozadaniowych liniowych układów regu- 
lacji automatycznej, opisanych w dziedzinie operatorowej se% kwadratowymi macierzami wy- 
miernych transmitancji, zapewniających jednocześnie: (i) kompletne dynamiczne rozsprzęgnięcie 
pomiędzy sygnałami zadanymi i wyjściami obiektu, (ii) uzyskanie z góry zadanych właściwości 
dynamicznych układu oraz (iii) eliminację na wyjściu w stanie ustalonym dodatkowych zakłóceń 
deterministycznych w obecności białych szumów gaussowskich. Rozważono trzy rodzaje struktur 
układu: z bezpośrednio mierzalnym wektorem stanu x obiektu oraz z wektorem stanu X od- 
twarzanym za pomocą stacjonarnego filtru Kalmana, przy czym odchyłki regulacji zdefiniowano 
odpowiednio jako e = yp—y i e=yo9—5. Wyniki badań zilustrowano wynikami symulacji dla 
układu nadążnego z wielowymiarowym liniowym obiektem ciągłym o dwóch wejściach i wyjściach. 
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CTAHHCJIAB BAHbKA, II93APbIl MOCKBA 


MCCJIEIOBAHHA MHOTOOYHKIIMOHAJIBHBIX CHCTEM 
YIPABJIEHMA IIPH IIOMONIM IIAKETA ASWUD 
NJIA AHAJIM3A M CMHHTEJA 
MHOTOMEPHDBIX CHCTEM YIIPABJIEHHA 


Pe3oMe 


B pa6oTe npe1CcTaBJIeHbl pe3yJIbTATbI HCCJIE1OBAHHA CJIOXKHbIX MHOTOQYHKHMOHAJIbHbIX CHC- 
TeM AaBTOMATHYUECKOTO YHpABJIEHHA, OIHCAHHbIX B OlepaTOpHOŃ OÓJIacTu se% KBAaHpaTHbIMU 
MaTpUIaMH pauNuoHaJbHbIX NepeAaTOUHbIX QyHKIIHŃ. CHCTEMBI yrpaBJIeHHA OJJHOBDeMEHHO OÓe- 
CIIEUHBAIOT: (i) KOMIIJIEKCHyFo NUHAaMHUECKYło pa3BA3KY MEXKJY 3AJNAHHbBIMH CHTHAJJAMM HM BbIXO- 
NaMM OÓbekTa, (ii) nojryueHie 3apaHee 3ANaHHbIX JIUHaMUUECKHX CBOŃCTB CHCTEM, a Take 
(iii) ycrpaHeHne Ha BbIXOJIE, B YCTAHOBUBIIIEMCA peXXUME, JOIIOJIHMTEJIbHbIX HETEDMUHHCTHYUECKHX 
TIOMeX IIDH HaJIMUMH OEJIbIX TAYCCOBCKHX IIIyMOB. PaCCMOTDpEHbl 3 BHNA CTPYKTyp CHCTEMBl: C He- 
NOCpeĄCTBeHHBIM H3MEpeHHEeM BEKTOPA COCTOAHHA X OÓbEKTA, A TAKXE C BEKTODOM COCTOAHHA 
£ BOCIIDOM3BOJHMMOTO C IIOMOIIIbIO CTANHOHApHOTO (bMJIEGTpA KaJbMaHa, NpH U€M OTKJIOHEHHE 
peryJlanMM OIIpej|eJJleHO COOCBETCTBEHHO KAK € = y,—y He =yq—y. HccJe1OoBaHHA HJIJIIOCTDPH- 
pyłoTca pe3yJIbTaTaMU MOJIEJIADOBAHHA JIA CJIEĄALIEŃ CHCTEMbl C MHOTOMEDHBIM JIMHeHHbIM 
HelIpepbIBHbIM OObeKTOM C JIBYMA BXOJJAMM M BbIXOJAMA. 
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Summary. A mapping of pseudo-Riemannian manifold (M, g) onto manifold 
(M, g) is said to be geodesic if it maps geodesics onto geodesics. A manifold (M, g) 
admits such mapping if and only if there exists a bilinear symmetric form a satisfying 
(3). In the paper it was proved that if an Einstein manifold admits such mapping, 
then its eigenvectors are non-isotropic. Making use of this result the local structure 
theorem was proved. 


Key words: Einstein manifolds — geodesic mappings. 


INTRODUCTION 


Assume that (M, g) and (M, g) are two n-dimensional pseudo-Riemannian 
manifolds of the class 6”. A mapping y: (M, g) > (M, g) is said to be geodesic 
if y maps geodesic lines of the manifold (M, g) onto geodesic lines of the 
manifold (M, g). Then their metrics are said to be geodesically corresponding. 

Suppose that both g and g are metrics on the same manifold M. Let 5(M) 
be the ring of differentiable functions and X(M) the 5$-module of differentiable 
vector fields on M. Each of the conditions below is necessary and sufficient for 
the metrics g and g to be geodesically corresponding: 


Y=FY+(XY)Y+(YyY)X, (1) 
Wx9)(Y, Z) =2(XV)g(Y, Z)+(YV)gl(X, Z)+(ZY)g(X, Y), (2) 
Wya)(Y, Z) = (Yp)g(X, Z)+(Zo)g(X, Y), (3) 


for all X, Y, ZeX(M), where o, V e $(M), V andV are the Levi-Civita connec- 
tions with respect to g and g, and a is a symmetric non-singular bilinear form 
on M. 
In a chart (U, x) on M the local components of g,g,a, Xp and XyV 
given by 
gu=g(X, X), du=J(X,X), aj=a(X,X), 


Pi=X;p, W,=Xy 
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satisfy 
d;; = exp(2Yy)9" gsi9ej» (4) 


j ) (5) 


where X, = 0/0x'eX(U) and the g” are the components of (g,,) *. 

If the geodesic mapping is not affine, then it is said to be non-trivial, which 
is equivalent to the condition Xy £0 (or Xo ź O). 

A manifold (M, g) is said to be an Einstein one if its Ricci tensor S$ is 
proportional to the metric tensor g, i.e. 


S(X, Y)=(n—1)Bg(X, Y), B=const, (6) 
holds for each X, YeXx(M). 


Ń i 
(B=, —exp(2Y)9" gsiV;» W; z 2(n+ 1) X; (ts 


PRELIMINARIES 


On the account of (1), the curvature tensors and Ricci tensors of (M, g) and 
(M, g) are related by 


R(X, Y)Z =R(X, Y)Z+P(X, Z)Y—P(Y, Z)X, (7) 
S(X, Y) =S(X, Y)+(n—1)P(X, Y), (8) 

where 
P(X, Y)=H;,(X, Y)-(XW)(TY) (9) 


and H, is a Hessian of the function Ye %$(M). It is well-known that the Weyl 
projective curvature tensor 


W(X,Y)Z=R(X, PZ-——[S(X, Z)Y-S(Y, Z)X] 


is invariant under geodesic mappings, ie. W = W. 
It is interesting that the class of Einstein manifolds is closed with respect to 
the geodesic mappings. We have: 


Theorem 1 [5] 


If the Einstein manifold admits non-trivial geodesic mapping onto a mani- 
fold (M, g), n> 2, then (M, g) is also Einstein manifold. 

Einstein manifolds of dimension < 4 which are not of constant sectional 
curvature do not admit non-trivial geodesic mappings. 

Let (M, g) be a Riemannian or pseudo-Riemannian manifold admitting 
non-trivial geodesic mapping onto a manifold (M, g). Then there exists 
a l-form do and a bilinear form a such that the relation (3) holds on M. 
We shall assume that the condition 


H,(X, Y)= Ba(X, Y)+ng(X,Y), X,YTex(M), 
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is (locally) fulfilled on M, where Hog is a Hessian of o and u, Be$(M). The 
condition (11) is equivalent to 


P(X, Y) = Bg(X, Y)—Bg(X, Y), (12) 
where P is as in (9) and Bet(M) ([6]). From (6) and (8) it follows that on 
Einstein manifolds admitting geodesic mappings the condition (12) is satisfied 
and B= zag e where r is a scalar curvature. 

n(n— 1) 
Proposition 1 [6] 

If Bź40 at a point peM then a vector © defined by g(X, $)=Xo, 
XeT,(M), is non-isotropic (non-null). 
Proposition 2 

If B40 at a point peM then a vector Y' defined by g(X, F)=XVv, 


XeT,(M), is non-isotropic. 


Proof 


Choose a local coordinate system (U, x) so that peU. Assuming that the 
vector Y' is isotropic (null), we have y,y' = 0, where y'* are components of the 
vector Y. By covariant differentiation, in virtue of (12), we obtain 


By, ck By'gii 


whence By' = By,g'". Hence, taking into account (4) and (5), we infer that g, 
and Y, are linearly dependent and the thesis follows from Proposition 1. 


EINSTEIN MANIFOLDS ADMITTING GEODESIC MAPPINGS 


Let a be a differentiable symmetric bilinear form on U, s M satisfying (3) 


1 t 

and having t different eigenvalues 4, ..., A. From very definition, at each point 
peU, they coincide with the eigenvalues of the eadomorphism 4, of tangent 
space T,(M) corresponding to a, ie. g(AX, Y) = a(X, Y) for all X, YeX(U,). 
Let (U, x) be a chart on M such that U £ U,. Suppose that b is an eigenvector 
of the matrix a;, corresponding to the eigenvalue A « =1,..., t, Le. satisfying 
the condition 

(a,; — Ag,)G = 0. (13) 


? ś i r ł = 
We define generalized eigenvectors u corresponding to A as follows 


a 1; a 
(a; — Ag;;)u* = Di, 


rl 


(a; — Ag.) = ZZ BAH (14) 


3 — Szczecińskie Roczniki Naukowe t. IX, z. I 
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Differentiating (13) and (14) covariantly with respect to x* and making use of 
(3) we have 


0d, + pig — AV, +(a;;— Agi) VA, =0, 
1 l; a 1 a 1; [74 
PU PU ik — Ax U; + (Qi; — Agi) Uk = 0z;k» 
r r: LJ r a r. r_1 
PU, + PW gk — AU; +(a,j—Agijuk= U zk (15) 


(the semicolon denotes covariant differentiation on (M, g)). From the above 
relations one can obtain ([2]) 


a; l 1;l-1 2;1-2 1-2; 2 
PpiVuy+Q,U u, OU w +...trQ,U Ut 


ja: - © l+la A; 
+pj wi, Poi h A uji =(. (16) 


If a manifold (M,g) admits non-trivial geodesic mapping determined 
by A z< constant, then, in virtue of (16), we have Di, z 0. On the other hand, the 
relation (14) results in U'U, =( for r=|1,...,l—1. Transvecting the second 
equation from among 1+ 1 relations (15) with "u' and the I-th equation with 5;, 


like (16), we obtain 
QV u ŁPU u PU u, +...+Q, U Ut 


1-1; 


a la - ©; 
+9, U zły UB! = 0. (17) 


1-1 


i 


© 


. a r . . ... —1d 
Hence, since vectors v, u are linearly independent and condition u, v 


: a 
holds, we have g,v' = 0 and g,u' = 0 forr = 1,..., 1-1. Because 2p= Y TA, 
=1 


where z, denotes algebraic multiplicity of 2 and 4 = A(x"=*"e), n, =0, 
Ng = Ty +Tą+...+Tp-1, E=2,...,t (L1, 2]), therefore from (16) we get 
6.=Ag, b=0 for jzi., (18) 
where i, =n,+1,..., h,+T,, AE F(U). 
Now we shall prove 


Theorem 2 


On an Einstein manifold (M, g) with B ź O all eigenvectors of the matrix 
a,;(p), peU,, are non-isotropic (non-null). 
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Proof 
Suppose that on the Einstein manifold with B ;+*0 the eigenvectors v 


a 
corresponding to the eigenvalue 4 of the matrix a;,(p) are isotropic. Differen- 
tiating (18) covariantly with respect to x* then transvecting the resulting 


equation with 0* and applying the relation Qi, 0 if jćj,, we obtain 


D,;xU' = (0. Therefore, from (11) and (13), we have 


Bi+u=0. (19) 


Moreover differentiating (5) covariantly with respect to x*, in view of (3), 
(12) and (11), we have 


u= —Bexp2y)+o'y,, 9'=o;g". 
Differentiating the above equation once more we easily get a, = 2Bg,, whence 
t 


u=B(); 1,4)+ C, where C = constant. Consequently, from (19) it follows 
a=l 

either B = 0 (and C€ =0) or the manifold (M, g) admits trivial geodesic map- 

ping which contradicts the assumption. This completes the proof. 


From the above Theorem it follows 


Proposition 3 


If the Einstein manifold with non-zero scalar curvature admits geodesic 
mapping, then the matrix a,;(p), peU,, of the form a is nondefective, 1.e. for 
each eigenvalue its geometric multiplicity is equal to the algebraic one. 


LOCAL STRUCTURE THEOREM 


Assume that on a manifold (M, g) the form a is nondefective and suppose 
that it has k different eigenvalues with multiplicity equal to 1 and t—k different 
eigenvalues with multiplicity greater than one. In the case under consideration, 
from papers [1, 2] and results of the previous section, it follows 


Proposition 4 


Let a be a form satisfying the condition (3) and suppose that for each 
peU, SM its matrix a;;(p) is nondefective, then in some neighbourhood of 
p there exists a coordinate system such that the metric tensors g and g take the 
form 


k t t t 
g9= X e.[1(6-407(x93+ X TI 4—/0"9:i, deedwie, 
u=1 bz o=k+1 fa 
u 
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ś= [1 "LZ 6607" [LUG/0"(bey+ 
"i pzu 


t 


+5 ("TI U,-f)"d:,, deedzet , (20) 
B=1 


o0=k+1 


Bu 
where J.= Je"), J, = const 40; e,= +1; | a= hę pW CL OIMESENE A 
M ZL Fq= l, 1, > Bi 4,,j, =h bl, RORŻy PAR PTY! nę =0, N, = Ty 
FT Fo FRA ="2 3,4016 aid dry”, ..., x"e"Toj are metric tensors 


2 
on 1,-dimensional submanifolds V. 


From (20) follows that M is locally warped reducible. Assume now that 
M is locally warped product manifold. Let pe M and U, be a neighbourhood 
of p such that U, = V,x,, V;x...x,, V„, 0, E8(%) and Vox, R*"x...x,_R' is 
the space of constant curvature K. Such decomposition of U, is said to be 
a K-decomposition whereas the manifold (M, g) is said to be of type V(K) 
([10]). 

As it is known, a manifold (M, g), n > 3, is of constant sectional curvature if 
and only if the Weyl projective curvature tensor W vanishes identically on M. 
A manifold (M,g) is an Einstein one if and only if g(W(X, Y)U. V) 
+g(W(X, Y)V, U) =0. If an Einstein manifold (M, g) admits geodesic map- 
ping onto (M, g), then the condition 


g(W(X, VU, V)+g(W(X, YJV,U)=0. X,Y,U,Vex(M), (21) 


is satisfied. 

Considering (21) in such coordinate system in which metrics g and g take 
forms (20), in virtue of [6, 10] and the above made considerations, we easily 
obtain 


Theorem 3 


An Einstein manifold (M, g) admits geodesic mapping onto an Einstein 
manifold (M, g) if and only if: 

(i) the manifolds (M, g) and (M, J) are of types V(K), K = —B,and V(K), 
K = —B, respectively, 

(1) for each peM there exist neighbourhoods U, = x, 4x... x, V» 
and U, = Vox, V;x...xg,V,„., 0, € %(V,), such that V, and W, are geodesically 


corresponding k-dimensional spaces of constant curvatures K and K, 
(ii) (V,, 9),0=1,2,..., m, are n,-dimensional Einstein spaces with scalar 
curvature r, = n,(n, —1)K,, where 


o” 


= = = ab 
=. Ję,a 0e,b 1 
40, 


sab 
_  Oga0eŁJ 


p: 40, 


K +Bo, = + BO, = const, 
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Jęa = Ka(0,),94” = (dw) ',a,b= 1,..., k, 4 and g being metric tensors on V, 
and V, respectively, 
(iv) GO Gd, zę 9%6,4. 


If an Einstein manifold (M, g) admits geodesic mapping onto an Einstein 
manifold (M, g) then the functions f, from (20) must satisfy the system (12). 
So, in virtue of Theorem 3 and the results of [10] (pp. 156—162), by an 
elementary calculation we obtain 


Theorem 4 


If an Einstein manifold (M,g) with $ =(n—1)Bg, B ;< 0, admits geodesic 
mapping onto an Einstein manifold (M, g) with S =(n—1)Bg and g, 40 at 
a point peM, then in some neighbourhood of p there exists a coordinate 
system such that the metric forms g and g take one of the following form 


l 
= zp = X) h.. dyxy'd i 
j EEG OR POSH 


MMNECH 


GAJ TH Oj, de dz, (22) 


l 
x (cy)" 
where c;, c, =const, h=h(x”,..., x") is the (n—1)-dimensional Einstein 


1 = 
metric such that S;,,, =(n—2)(c,B—c.)h,,,, B=(c,)"ca, I4, ją = 2,..., n, 


l 1y2 1 ji 4. 
jg = ke xa = a + (c, —x')(c; —c,)? CEE dx 
+ (c, — x')(c, —c,)" A dx'?dx72, 
Jj= g (5 OESIE x 
EARZAECIA 5. (Eg "> (eg) 


X gi, jj dx'"dx'' + dx?dx*2, (23) 


z (ee(ee 1a 


1 1 
where c;, c, = const, h=h(x*,...,x''*') is the z,-dimensional Einstein 


1 h 2 2 
metrie such that S;,,, = (c, —DB(Ev=c.)"" hy, > AOF"*, ....m) is-the 


tą-dimensional Einstein metric such that SA (1, —1)B(c,—c,)"* lisy. 


laj 2 awe PID JT, FA 4a ly Pr, SKG B= Ble)" la)?" >, 


("x") 


t k 
(dx)? + 2 T[(/,—x)A,J;,,dxedxie, 


u=1 LAN: Q(c*) o=ktlu=l 
ki 1 t k Key] k (x — x") z 
z = x AL (0-7 p » II Q(x") 8 du 


v£u 

t k 
F2 (ŻUI (i-*)Abidxsash), (24) 

o=k+1 u=l1 


t 


where A,= [| (—/.)”, f,,f, = const40, o,w=k+1,...,t; Q(z) is 


o=k+l 
OE 
a polynomial of the form Q(z) = (—1)**'4Bz**'+Cz' + Az-1z* !'+...+A;z 
t 
+4C;,C,=B [| (%) % €, 4y,..., 4;- , = const, for each gf, are a root 
o=k+1 
of the polynormnial OQ(z),k > I,t <2k+1, J are metric tensors of z,-dimensional 
1 
m (Gz= DK, di, and K, =; 4(7077'4,0 (): 


o 
Einstein manifolds with Ś;;, 


sę > M tlywee RZEP NY FAN = Tr Tt ha plyny = 290960, 251. 
t 
The function Y is defined by exp(—2y)=x'...x* [[ (f,)*. 
o=k+1 
Remark 


In the case B = 0 the local structure of Einstein manifolds admitting geo- 
desic mappings can be easily obtained on the ground of [8]. 


PROJECTIVE MOTIONS IN EINSTEIN MANIFOLDS 


A vector field X on a manifold (M, g) is said to be a projective motion if the 
one-parameter local group of transformations generated by X in a neighbour- 
hood of each point peM leaves geodesics unchanged. In order that 
a one-parameter group of transformations generated by a vector field X E X(M) 
be a group of projective motions it is necessary and sufficient that X satisfies 


L;g=h, 
W„h)(Y, Z) = Żg(Y, ZVY +gl(V, ZJYY +g(Y, V)ZY, 
Y,Z,Vex(M), Ve5(M), 


where Ł, is the Lie derivative in direction X. 
A projective motion is said to be an affine motion if Y = const. We have the 
following 
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Theorem 5 [5] 

If an Einstein manifold with non-zero scalar curvature admits non-trivial 
geodesic mapping, then it also admits non-trivial projective motion. The con- 
verse of this theorem is also true. 

If an Einstein manifold with vanishing scalar curvature admits non-trivial 
geodesic mapping, theń it admits an affine motion. 


STANISŁAW FORMELLA, JOSEF MIKESH 


GEODEZYJNE ODWZOROWANIA ROZMAITOŚCI EINSTEINA 


Streszczenie 


Geodezyjnym odwzorowaniem rozmaitości pseudoriemannowskiej (M,g) na rozmaitość 
(M, g) nazywamy odwzorowanie, które zachowuje linie geodezyjne. Warunkiem koniecznym i do- 
statecznym na to, aby rozmaitość (M, g) dopuszczała takie odwzorowanie, jest istnienie na niej 
dwuliniowej formy symetrycznej a spełniającej warunek (3). W pracy wykazano, że na rozmaitości 
Einsteina dopuszczającej geodezyjne odwzorowanie wszystkie wektory własne formy a są nieizo- 
tropowe i na tej podstawie udowodniono twierdzenie o lokalnej strukturze takich rozmaitości. 


CTAHHCJIAB WOPM3JIJIA, E.CJO MMUK3JULI 


FEOTE3AUECKOE OTOBPAXEHME 
MHOTOOBPAJMA JDAHIITEAHA 


Pe3roMe 


Feoz1e3uueckuM OTOOpaxeHieM MHOrooOpa3na (M, g) Ha (M, g) Ha3biBaeTca oTo6paxeHie, 
KOTOpoe CcoXpaHaeT TeO|Ee3HUECKHE JIMUHHH. UTOÓbI MHOrooOpazne (M, g) NorryckaJio TaKOE OTOÓ- 
paxeHHe, H£OÓXOJIAM5BIM H NOCTATOUHbIM YCJIOBHEM ABJIAETCA CYLIIECTBOBAHHE Ha He NBYJIMHeh- 
Hof cuMMeTpHueckoń (DOpMBl a YHOBJIETBOPAEOLNEŃ YCJIOBHŁO (3). B pa6oTe npe4CTaBJIeHO, UTO Ha 
MHOTOOOpa3uu DKHIITeAHa, NONYCKAIOLNEE TEOJJE3HUECKOE OTOÓPAXEHHE, BCE COÓCTBEHHbIe BEK- 
TODpbI POPpMbI a ABJIAKOTCA HEH3OTDPOIIHBIMH. Ha 3TOlŃ OCHOBe NTOKAa3AHa TEOpeMA O JIOKAJIŁHOŃ 
CTpykType TakHuX MHOTOOOpa3NA. 
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Summary. The eddy-current method of nondestructive testing (NDT) involves 
the introduction of eddy-currents in a conductive test object by a time-varying field 
and the detection of the resultant field. The theoretical analysis of the electromag- 
netic field distribution in eddy-current systems is very difficult and unsolved at this 
time. In this paper the 2-D and 3-D computer simulation for eddy-current nonde- 
structive testing of metal objects is described. The results presented show that the 
finite element method can be successfully used for modelling field/flaw interactions 
and for providing training data for automated signal processing equipment. The 
results of calculations are expressed in terms of impedance diagrams and signal 
trajectories for various kinds of eddy-currents probes. 


Key words: finite element method — nondestructive testing of materials. 


INTRODUCTION 


The eddy-current method is widely used for problems of material science, 
which are, generally speaking, those of 

— measurement of electric conductivity at higher temperatures without 
using electrodes, 

— determination of the action of external and intercrystalline corrosion 
with a sensitivity of 10 ©, 

— gquantitative appraisal of the variation of electric conductivity and the 
cross-sectional area of the sample for plastic and elastic strain deformation, 

— investigation of the influence of hardening and ageing on the electric 
conductivity and the diameter of the sample, 

— recording first stress cracks and measurement of the increase in depth of 
the cracks, 

— occurrence of irreversible variations in conductivity and diameter of the 
sample during a fatigue test, 

— observation of diffusion phenomena connected with multi-layer pro- 
blems. 
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Tests of plates, shells and tubes are particularly important, in view of the 
increasing demand for materials subjected to high temperatures. They are used, 
for instance, for aerodynamic surfaces of space vehicles, power turbines etc. 
Heat-resisting alloys with silicide coating can resist temperatures reaching 
1400”C. Defects of the coating may be catastrophic for the entire structure, 
therefore nondestructive tests of such a structure before and during operation is 
essential. 

Eddy-current nondestructive testing systems detect defects in a conducting 
material by using a sensor which induces currents in the material and then 
observing the impedance changes of the sensor in the vicinity of the defect. In 
this paper the finite element method was applied in the development of a model 
for impedance calculations of probes in nondestructive testing. Because the 
dimensions of the flaws which are to be detected are small and the frequency of 
excitation is relatively high, one has to take a very fine adaptive mesh in 
a region near the defect to obtain an accurate solution. 


FINITE ELEMENTS SIMULATION IN NONDESTRUCTIVE TESTING 
OF MATERIALS 


The aim of using numerical models is 

— to describe the interaction between the variable excitation field, the 
induced eddy-current and the material flaws, 

— to furnish examples difficult to be reproduced by experimental means, 

— to obtain diagrams of the output quantities such as the difference im- 
pedance for various defects and geometrical forms of the object without the 
necessity of expensive experimental studies, 

— to serve as an assistance in designing and optimizing of probes for 
particular purposes. 

Numerical methods, those of finite elements in particular, offer some new 
tools for the analysis of eddy-current flaw detection systems. The problem of 
nonlinearity of the medium and its geometrical form has been solved in the 
case of two-dimensional and axially symmetric systems, the problem of analysis 
of three-dimensional systems remaining still unsolved, however. In addition to 
the difficulties connected with uniqueness of solution or the way of determining 
the boundary conditions there is a serious obstacle of limited capacity of the 
memory and calculation speed of computers. 


THE 2-D MODEL OF THE SENSOR FOR SURFACE TESTS 


The Fórster flaw detector is equipped with a number of probes for testing 
surfaces. They differ by their dimensions, and the selection of a suitable probe 
depends on the expected size of the flaw and the quality of the surface to be 
tested. A sectional view of the model of the Fórster probe is presented in Fig. 1. 
The model is axially symmetric, this fact being made use of only for determining 
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the frequency characteristics of the probe in the case of no flaw being present, 
however. The probe is composed of a screened coil with a ferromagnetic core. 
It has been placed on an electrically conducting plate containing a flaw to be 
detected. In order to reduce the error due to the analyzed region truncation, 
infinite elements have been introduced on the boundary. 


"2 Infinite elements 


ze and screen 
=1000, 9=0 


infinite 
elements 


conducting plate 
4=10 Ś/m, H=Lo 


Fig. 1. The model of the Fórster probe 


The calculations were performed by the method of finite elements, using the 
magnetic vector potential A defined as B = rot A. Similarly to the currents 
generating the magnetic field, this vector potential has only one component, 
A„. The distribution of the potential is described by the equation 


0 ( (4400) żę O(A,r) 


AL A 2 = | 
Or ru Oz * (U) 


ru  Ór Oz 


for the air, the core and the screen, 


0/1 6) 0 PAD) M, 

(a Or +a(x 0z ueiże (2) 
0 41 0(45r) [> 41 0(45r) CF 

A Ór +z(; 2) jon, (3) 


for the conducting plate. 


for the coil and, 
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Fig. 2. The element mesh used for analysis 


Figure 2 shows the finite element mesh used for analysis. This element mesh 
was obtained by means of a generator on the basis of macroelements. Next the 
mesh was subjected to an adaptive process which resulted in considerable mesh 
refinement inside the conducting plate in the neighbourhood of the flaw. The 
mesh generators and adaptive algorithms are part of the SONMAP-micro 
software package and have been described in [1]. 
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Fig. 3. Frequency characteristics of the probe of the Fórster defectometer: a) a mesh with 311 
nodes, b) a mesh with 587 nodes, c) adaptive mesh with 497 nodes 
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The H 2.835 Fórster defectometer analyses only the real part of the probe 
impedance. To determine the optimum operation frequency of the probe, the 
frequency characteristics of the probe resistance variation were determined. 

The left-hand figure illustrates the relative variation of the resistance of the 
probe during the detection for a surface flaw of a depth of 0.5 mm and a width 
of I mm. The detection sensitivity for such a flaw increases with frequency. It is 
seen that the results obtained for the fine mesh and the adaptive mesh are, in 
this case, similar. The frequency characteristics for a flaw detecting below the 
plate surface is represented in the right-hand figure. The dimensions of the flaw, 
which was located at a depth of I mm below the surface, were 1 x 1 mm. Such 
a defect cannot be detected, if the feed frequency of the probe is too high, the 
advantages of the use of the adaptive method are therefore obvious. 

The characteristics presented above have been obtained for a location of 
the flaw corresponding to the output signal of maximum level. The resistance 
variations of the probe are obtained during the passage of the probe above the 
flaw. Their form shows the size and the character of the flaw. Figure 4 shows 
such characteristics, obtained by simułation, for two frequencies and different 
element meshes. The defect was located at the surface and its dimensions were 
not changed. 
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Fig. 4. Relative resistance variation of the probe during passage above the flaw: a) a mesh of 311 


nodes, b) a mesh of 587 nodes, c) an adaptive mesh of 497 nodes 


The only fact considered in the analysis presented here is the variation of 


the probe resistance, because the defectometer described here analyses the 
variation of only the real part of the impedance of the sensor. It has been 
shown by computation that the induction reactance dominates the resistance, 
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its variation under influence of the presence of the flaw being insignificant, 
however. This is due to the fact that the structure of the sensor is such that 
most of the field lines are closed inside the probe so that they do not pass 
through the material (Fig. 5). The algorithms described enable us to optimize 
the structure of the sensors. 
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Fig. 5. Field lines of the sensor 


AXISYMMETRIC MODEL FOR TUBE TESTING SIMULATION 


The model under investigation is shown in the figure below (Fig. 6). It 
consists of the differential coil inside a tube. The material defect is also shown. 
This model is adequate to the test sample, which is described by DIN 54141. 
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Fig. 6. The differential coils inside the tube with a flaw 
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The use of the cylindrical two-dimensional models can be a good 
approximation of three-dimensional problems. Similarly to the model of 
the sensor for surface tests the calculations were provided for the mag- 
netic vector potential A. The feed frequency is low enough and the displace- 
ment current term in Maxwell's equations is negligible. For sinusoidal 
excitation the eddy-current field equations become the form (1), (2) 
and (3). 

The flux density lines nearby sensor obtained by FEM are shown in the 
Fig. 7. The feed frequency is equal to 50 kHz. 


Fig. 7. Flux density lines 


For purposes of eddy-current testing, the most important feature of 
the FEM analysis is the ability to compute the impedance at the ter- 
minals of the sensor with respect to defects. A plot of Z on the 
complex impedance plane for various positions of the probe as it moves 
past a defect gives the predicted eddy-current probe signal trajectory for 
that defect. The impedance trajectory can serve as a certificate of the 
product. 

The exemplary impedance trajectories for outside and inside defects are 
shown in Fig. 8. 

The problem of choosing a proper value of frequency is of great import- 
ance. The optimal frequency of the input signals depends on the tube wall 
thickness and defect depth. The output signals for various frequencies and for 
arbitrary position of the sensor with respect to the inner defect are shown 
in the Fig. 9. 


48 


outer flaw (x 0.0008) inner flaw (x 0.013) 
1.07 AX 


Fig. 9. Normalized impedance diagram 


THE THREE-DIMENSIONAL NDT MODEL 


The structure of a probe of an eddy-current defectometer being axially 
symmetric, it can easily be described by means of a two-dimensional model. 
Considerable difficulties arise, however, if material flaws, which are not sym- 
metric, in general, are to be modelled. In this case a three-dimensional model 
must be used. In addition to the difficulties of mathematical formulation of the 
problem there are also difficulties of numerical nature, due to the computing 
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power of the existing computers, being too low in the case of analysis of vector 
fields, in particular. Below, we shall present an idea of determining 3-D fields 
by means of a scalar function. The analysis will be made by using the magnetic 
scalar potential $ and the electric vector potential T. They are related with the 
magnetic field intensity vector H and the current density J by the equations 


H = T-gradó and J=rotT. 


The distribution of the potential Tis a known function describing the 
excitation current density. This method excludes the existence, in the region 
analyzed, of unknown eddy-currents, which cannot be described under such 
conditions by means of the potential T. Since the principle of measurement with 
Fórster defectometer is based on the interaction of eddy-currents with the 
probe, special impedance boundary conditions may be used to take them into 
account. The model analyzed is of a relatively simple configuration because it 
constitutes a rectangular coil located above the conducting plate, Fig. 10. 
A cross-section of this model in the x, y-plane is shown in Fig. 11. 


Fig. 10. Rectangular coil above a plate with a flaw 
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Fig. 11. The x, y-section of the model 
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A three-dimensional element mesh was obtained by means of layer-gene- 
rator [1]. It extends the existing triangular mesh in the direction of the 
third dimension, forming first prisms, then tetrahedrons. The application 
of tetrahedrons as finite elements has many advantages. The shape functions 
are simple and integrable analytically. It makes possible defining of stiffness 
matrix coefficients explicit. Next, such a mesh can be refined by means of an 
adaptive method, as a result of which new tetrahedrons are formed [1]. A draw- 
back of the use of tetrahedron elements is that a considerable number 
of such elements are needed for a mesh. The initial mesh had 405 nodes and 
1536 elements. 

If the scalar potential $ is used the currents exciting the magnetic field are 
taking into account by considering the electric potential T. The current density 
in the coil has only two components, J, and J,, therefore there is only one 
component 7, (Fig. 12). The potential T exists inside coil only. 


A 1 [Ałem 
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Fig. 12. The value of T,(x, y) inside the coil 


The conducting plate below the coil was taken into account by imposing an 
impedance condition at the boundary of the region 


0% 3ó 06 u,(1 —j)ó 2 
——=x|—+— here: = ——— = |-— 
dy (38 SJ OE 4 o 3 AŻ 
U, — relative permeability of the medium. 


Taking into account the impedance boundary condition results in the 
additional terms that are to be introduced into the stiffness matrix of the 
element. 

The model described here has many symmetry planes so that only 
one-eighth of the region can be taken for computation. The zero Dirichlet 
condition was imposed at the z = O plane and the zero Neumann condition in 
the remaining planes of symmetry. 
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Fig. 13. The infinite element 


Because a truncation of the region may be a cause of considerable errors, 
infinite elements were used [2]. Such an element, made to suit the tetrahedrons, 
is shown in Fig. 13. If the scalar potential of the magnetic field © in an infinite 
element is to decrease in agreement with the imposed function / (r), it suffices to 
assume for $ the following expression 


3 
$'(r, O, p) a A > N;Q;, 
i=1 
(z zz) 


with the classical shape functions for triangular elements, N,. For the function 
f (r) we assume, in the simplest case, / (r) = 1/r”, n > 1. Then, the coefficients of 
the matrix can be determined by analytical means. 
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Fig. 14. Section of 3-D finite element mesh for z =0.125 cm 
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Fig. 15. Section of 3-D finite element mesh for x = 0.06 cm 


Figures 14 and 15 present sections of element meshes which were 
used for computation (1400 nodes), showing mesh refinement obtained 
by adaptive methods. Those meshes were used for simulating the detection 
of a surface flaw. The flaw was modelled by a rectangle of the dimensions — 
0.125x0.125 cm. It was assumed that the magnetic field is perpendicular 
to the plate surface at the rectangle modelling the flaw. This assumption 
was realized by means of the impedance boundary condition with the parame- 
ters u, = 1000 and y= 1 S/m. The flaw was translated along the x-axis 
of the probe by steps of 0.125 cm, the number of which was eight. Relative 
variations of the coil resistance for different locations of the flaw are shown 
in Figs. 16 and 17. 


Fig. 16. Relative variation of coil resistance for three meshes, f = 1 kHz: (1) initial mesh, (2) fine 
mesh, (3) adaptive mesh 
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Fig. 17. Relative variation of coil resistance for three meshes, f = 100 kHz: (1) initial mesh, (2) fine 
mesh, (3) adaptive mesh 


The simulation process was carried out with two frequencies, namely 
those of I kHz and 100kHz. It should be stressed that the impedance 
boundary condition works correctly by small penetration depth of the 
field, which means high frequencies. The tests which were performed 
showed that there was agreement within the limits of 10% with other methods 
of power losses analysis in the conducting material by solving the Helmholtz 
equation. 


CONCLUSIONS 


A number of methods of computer simulation of eddy-current flaw 
detection systems have been described. As regards 2-D cases, models 
of axial and plane-parallel symmetry have been discussed. Solutions 
were based on the Helmholtz equation for the magnetic vector poten- 
tial A. 

The three-dimensional model is based on the scalar potential q. 
An essential novelty is the use of an impedance boundary condition 
for modelling eddy-currents in the conducting plate and for modelling 
a flaw. The use of infinite elements for tetrahedrons, first published 
in [2], reduces considerably the error caused by cutting-off part of the 
region. 

The examples shown prove the usefulness of numerical simulation of 
eddy-current flaw detection systems for the design and optimization of sensors. 
They are also applied in the domain of modern microcomputer aided defec- 
tometers. 
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KONSTANTY M. GAWRYLCZYK 


ZASTOSOWANIE METODY ELEMENTÓW SKOŃCZONYCH 
DO SYMULACJI UKŁADÓW DEFEKTOSKOPII WIROPRĄDOWEJ 


Streszczenie 


Metoda niedestrukcyjnego badania wiroprądowego polega na wprowadzeniu do przewodzą- 
cego obiektu prądów wirowych za pomocą zmiennego w czasie pola elektromagnetycznego, a na- 
stępnie detekcji powstałego pola. Analiza teoretyczna rozkładu pola w systemach wiroprądowych 
jest trudna i dotąd nierozwiązana w sposób kompletny. W pracy tej opisano dwu- i trójwymiarowe 
komputerowe symulacje badań wiroprądowych obiektów metalowych. Otrzymane rezultaty wyka- 
zują, że metoda elementów skończonych może być skutecznie użyta do modelowania oddziaływa- 
nia pola z defektem oraz do tworzenia bazy danych dla aparatury rozpoznającej sygnały od 
defektów. Wyniki obliczeń przedstawiono w postaci wykresów impedancji oraz trajektorii sygnału 
dla różnych typów sond wiroprądowych. 


KOHCTAHTBI M. FABPBIJIbUbIK 


IIPAMEHEHME METOJTA KOHEUHbBIX EJJEMEHTOB 
AJIA HMHTATOPA 
CHCTEM BMXPEBOTOKOBOMŃ JIEPEKTOCKOIIMMA 


Pe3roMe 


MeToq1 Oe31eCcTpyKTUBHOTO BHXPpEBOTOKOBOTO HCCJIETOBAHHA COCTOHT B TOM, UTOÓBI BBECTH 
B HIpOBOĄAINMA OÓBEKT BUXpEBbBIE TOKH C IIOMOILNbIo IIEpeMEHHOTO IIO BDEMEHH MATHHTHOTO NOJIA, 
a IOTOM J1eTeKTHPpOBaTb OOÓpa3OBABIIeeCcA IOJIe. TeopeTuueckHui aHaJIH3 pacHperejjeHHA NOJIA 
B BHXpEBOTOKOBBIX CHCTEeMaX OUeHb TDYHHBIŃ H AO CHX IIOp He peLIeH HeJlAKOM. B 3Toń pa6oTe 
ONMCAH ĄBy- M TpexMEpHkIŃ KOMIIBIOTEpHbIi MMHUTATOp MIA BHXPEBOTOKOBBIX HCCJIETOBAHH 
METAJIJIMUECKHX OÓBEKTOB. I[IOJTyueHHbie pe3yJIbTATbI IOKA3BIBAŁOT, UTO METOJ KOHEUHbIX 3JIEMEH- 
TOB MO%XEeT 3QheKTUBHO IIDAMEHATBCA JIJIA MOJIEJJADOBAHHA BO3JIEŃCTBHA IIOJIA C Te(DEKTOM M JIA 
CO3NAHHA Óa3 NAaHHBIX NIA ANIApAaTypbI, pACHO3HAIOINEH CHTHAJIBI OT „AeQeKTOB. Pe3YJIbTATbl 
pacueTOB IIpe|CTABJIEHbI B BHJe TpaQHKOB UMIIEJEHCA HM TpAaeKTODHH CMTHAJIA NIA PAZJIMUHBIX 
BHNOB BHXpEBOTOKOBBIX 30HJI. 
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Summary. The paper reviews possible applications of static magnetic and 
Electron Paramagnetic Resonance measurements in studies of high temperature su- 
perconductors and related compounds. Particular examples of these applications are 
given and discussed. It is shown that magnetic measurements reveal the existence of 
different magnetic phase transitions, whereas the EPR measurements yield unique 
information about the coupling of a variety of paramagnetic centers with their crys- 
talline environment. 


Key words: superconductors — static magnetic measurements — Electron 
Paramagnetic Resonance — transition ions — electric conductivity. 


INTRODUCTION 


Since the discovery of high temperature superconductivity (HTS) in cera- 
mics [2, 18] great effort has been done to understand the origin of this phono- 
menon. A great deal has been learned about magnetic properties, crystal struc- 
ture, phase composition and the electronic structure of magnetic ions being 
present in superconducting and related compounds. Many experimental tech- 
niques were applied (some of them are very much sophisticated) in order to get 
the information useful in the microscopic interpretation of the HTS. 

In this paper we shall review some results of our experimental studies of 
magnetic properties of the compounds under discussion derived from static 
magnetic and EPR methods. We shall also try to interpret some mechanisms 
of conductivity. As deduced from the results of magnetic measurements 
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we were able to discover many unusual phase transitions occurring in phases 
related to superconducting ones. The EPR measurements supplied important 
information about the electronic structure of different magnetic ions and their 
interactions with the surrounding. 

The results described in the paper are reported in many of our previous 
publications and therefore it may be difficult, even for a specialist, to work out 
a general view about the conclusions resulting from them. This paper is intend- 
ed to help a reader interested in superconductors to overcome these difficul- 
ties. 


INFORMATION ABOUT EXPERIMENTAL METHODS 


AII ceramics studied by us were prepared by the standard solid state reac- 
tion technique. Appropriate amounts of Re,O; (Re stands here for a rare-earth 
element), Cu and BaCo; were carefully mixed and pelleted. The pellets were 
then annealed in flowing oxygen at the following four step reaction: 


920*C(20hs) > 925*C(20hs) > 930*C(30hs) > 940*C(20hs). 


After this procedure the samples were powdered and recalcinated at 
940”C(12hs). In order to increase the content of oxygen the pulverized samples 
were kapt in a hot furnace (920”C) with flowing oxygen. When thermal equili- 
brium was reached the temperature was reduced to 450'C (6hs) and then the 
samples were kept for other 6 hours at 400?C. In order to obtain the tetragonal 
phase, a part of the oxygenated powders was annealed at 650 C (6hs) in flowing 
He and then cooled fast to room temperature in a reducing atmosphere. 
The phase composition and structure characterization of the obtained 
samples were done using a Philips X-ray powder diffractometer with a Co-K, 
cathode. An example of the X-ray diffraction pattern (NdBa,Cu;O7 _5) 
for superconducting (orthorhombic) and non-superconducting (tetragonal) 
phases is shown in Fig. 1 [5]. For this particular case the major phase in 
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Fig. 1. X-ray diffraction pattern for orthorhombic (a) and tetragonal (b) phases of NdBa,Cu;O,_, 
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high temperature superconducting sample could be indexed on a base of an 
orthorhombic unit cell with the lattice constants: a = 0.38661 (5) nm, 
b =0.391632(2) nm and c= 1.17607(9) nm. 

The major phase in the non-superconducting sample could be in- 
dexed on a basis of a tetragonal phase with a=b=0.38915(4) nm and 
c=1.18390(10) nm. 

The magnetization measurements vs. temperature were carried out using 
a PAR 155 vibrating sample magnetometer in the temperature range from 
4.2 K to 300K. The applied magnetic field induction was usually 1 or 
2x107*T. The measurements were carried out both in zero-field-cooled (ZFC) 
and field-cooled (FC) modes. 

The EPR measurements were done using a Varian E-4 X-band spectometer 
working with a 100 kHz frequency modulation of the steady magnetic field. 
EPR spectra were taken in the first derivative mode. DPPH standard was used 
as a frequency marker and the magnetic field was scaled using a precision 
NMR magnetometer. 


STATIC MAGNETIC MEASUREMENTS 


For the samples in orthorhombic symmetry the magnetic measurements are 
usually performed in order to obtain quick information about the transition 
temperature to the superconducting state by the observation of the Meissner 
signal. 


magnetizat ion (e.m.u./g) 


o 20 40 60 80 100 
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Fig. 2. Temperature dependence of the magnetization for GdBa,Cu;O,_, 

An example of the magnetization vs. temperature for such samples is given. 
in Fig. 2 for a particular case of the GdBa,Cu;0;_; sample [6]. For this 
sample the 7, outset was 83 K. From the FC line and the ratio L= —4nM/H 
one can calculate the percentage of the superconducting phase estimated in this 
case to be about 13%. Such measurements are standard ones and therefore 
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need not be described here in detail. The results for the samples with tetragonal 
symmetry are much more interesting. This phase is non-superconducting and 
normally exhibits an antiferromagnetic ordering below T„. However, Masuda 
et al. [15] reported that for LaBa,Cu;0-_„ prepared in a reducing atmos- 
phere one detects, at low temperatures, an anomalous behaviour of the inverse 
susceptibility vs. temperature. The origin of this behaviour was axplained to be 
due to the existence of an impurity phase La,Ba,Cu,O;, (4:2:2 phase). How- 
ever, we found [7] that for YbBa,Cu;O, similar anomalies are detected al- 
though the data derived from XRD and EPR measurements excluded the 
existense of the 4:2:2 phase. 
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Fig. 3. The plot of the inverse magnetic susceptibility versus temperature for YbBa,O-0, 


The plot of the reciprocal static magnetic susceptibility for YbBa,Cu;O, as 
a function of temperature is shown in Fig. 3. Three main regions can be distin- 
guished from the plot. From high temperatures down to 60 K, x " obeys the 
Curie-Weiss law. The line extrapolated from these data crossses the abcissa at 
a negative temperature © = —25 K indicating the antiferromagnetism of the 
system. The effective number of Bohr magnetons determined from the slope is 
3.85. In the region between 60 and 30 K a clear deviation of x * from the 
linear run was detected. From 30 K down to 4.7 K one again detects a linear 
run of x " but the line extrapolated crosses the abcissa at about I K which 
indicates a positive divergency of the magnetic susceptibility at a finite tem- 
perature. This would mean that below 1 K the substance exhibits an ordering 
of the ferromagnetic type. This suggestion was partially supported by the re- 
sults obtained from EPR measurements of Cu** ions. It proved that below 
20 K one detects a gradual increase in the value of g. At 4.7 K this value 
reached as much as 2.205(3). At the same time when approaching the lowest 
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temperatures investigated from above we noticed a phenomenon of a substan- 
tial increase in the linewidth. All these facts can be interpreted within the frame 
work of the fluctuation theory for magnetization (see e.g. [16]). 

Other examples of possibly new magnetic phase transitions in tetragonal 
phases of ReBa,Cu;O05 +, (where Re = Nd, Sm or Tm) were given in [3]. This 
paper gives clear evidehce that for all investigated samples a pronounced 
change in the slope of x * versus T is observed. The plots exhibited a positive 
divergency which again indicated a ferromagnetic behaviour. 

Interesting data about magnetic properties of EuBa; _ „Eu,Cu-0;_; in tet- 
ragonal phase are reported in [8]. The plot of the reciprocal statistic magnetic 
susceptibility for this compound is presented in Fig. 4. 
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Fig. 4. The plot of the inverse magnetic susceptibility versus temperature for EuBa,_ Eu _Cu,O,_, 


As results from the Hund rules, the ground state of the 4f* configuration 
of Eu** is 'F, implying that the magnetic moment u should be zero. However, 
one measures the experimental value of u equal to 3.4u5p. This behaviour can 
be very well interpreted taking into account the fact that the first excited 
state "F, has the energy equal to only 4A/k, = 350 K above the ground 
state being thus quite considerably populated at higher temperatures [18]. For 
EuBa,Cu3O0>_; a paramagnetic susceptibility was reported by Xiao et al. [19]. 
They could, indeed, explain their data based on the above-mentioned fact. Our 
magnetic susceptibility measurements indicate a Curie-Weiss behaviour above 
200 K. The effective magnetic moment ur obtained by least-square fitting 
is 4.10u5 per formula unit. This value is clearly different from that expected 
for Eu? * ions including contributions from thermally excited states. The dis- 
crepancy can, however, be explained assuming that a small fraction of Ba** 
ions is substituted by Eu?* ones with per = 7.94 Mg. This explanation was 
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also supported by the EPR data. From the contribution of Eu? * to the total 
magnetic susceptibility x = 0.07 was estimated. 

Below 200 K the plot deviates from the Curie-Weiss law and below 50 K 
one detects that there exists a strong decrease in x " with the temperature 
decrease. The' behaviour below 50 K may be explained by the crystal field 
effect. 


EPR MEASUREMENTS 


EPR measurements are mainly devoted to the investigations of two pro- 
blems: the Josephson effect and the electronic structure and local symmetry 
of paramagnetic ions in HTS and related compounds [17]. Due to the intrinsic 
sensitivity of the method a great deal of important information enabling a bet- 
ter understanding of the nature of superconductivity can be derived. In this 
chapter we shall limit ourselves only to the discussion of chosen aspects of the 
EPR application in studies of HTS and related compounds paying particular 
attention to the electronic structure and local symmetry of paramagnetic ions 
involved. 

Let us begin with a short discussion on the results obtained for Cu** ions. 
Typical EPR spectra for these ions in orthorhombic and tetragonal phases are 
shown in Figs. 5—7 for the case of TmBa;Cu30;_;. The signal presented in 
Fig. 5 clearly results from Cu?" ions feeling the action of the crystal field of 
rhombic symmetry. The components of the g-tensor were estimated to be: 


g, = 2.099(5), g, = 2.029(5) and g, = 2.235(5). 


Below the transition temperature 7, = 91 K basically the same line was observ- 
ed and slight changes detected (see Fig. 6) could be ascribed to the increased 
conductivity of the sample. 
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Fig. 5 EPR spectrum of Cu** for TmBa,Cu;O,_, in orthorhombic phase taken at room 
temperature 
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Fig. 6. EPR spectrum of Cu** for TmBa,Cu,O._, in orthorhombic phase taken at 77 K 
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Fig. 7. EPR spectrum of TmBa,Cu,O,_, in tetragonal phase 


Similar pattern was seen for many other samples indicating that Cu?* 
ions detected by EPR are not directly involved in superconductivity Fig. 7 
presents the EPR spectrum for the tetragonal phase of the compound. As seen 
from the figure the signal reflects an axial local symmetry. Thus, we may 
claim that although Cu** ions observed by EPR are not involved in supercon- 
ductivity nontheless the EPR spectra of these ions identify unequivocally 
the type of the sample ie. they differentiate superconducting and non-super- 
conducting samples. 

Temperature dependences of the EPR spectra of Cu** ions suppied also 
additional information about the exchange interaction between these ions. 
As seen from Figs. 8 and 9 at low temperatures, typically below 25K, 
one observes a relatively symmetrical line, the intensity of which increases 
strongly with the temperature decrease. It is suggested that this line may 
result from the pair interaction, involving oxygen bridges, of divalent 
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Fig. 8. Low temperature EPR signals due to exchange coupled pairs of Cu?* for 
Pr, sDy,.„Ba;Cu;O;,_, 
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Fig. 9. Low temperature EPR signals due to exchange coupled pairs of Cu?* for 
Gd, „DyBa,Cu.O0,_; 
copper — copper ions. [his interaction can be described by the following 
spin-Hamiltonian [1]: 
X =g,U5B:S,+galgB':S,+JS4' 83, (1) 


where J is the exchange integral and the rest of the symbols have their usual 
meaning. 
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If one assumes an isotropic exchange interaction then the energy W of the 
triplet state is: 


W =J/4+1/2(g, +9»)u5B (2) 


and it is different from that corresponding to the singlet S= 0 state. For 
example, for the case of GdBa,Cu30-_; a single EPR line (shown in Fig. 8) 
with the g-factor equal to g = 1/2(g, +g>) = 2.27 was observed [9]. 

The value of J can be estimated from the temperature dependence of the 
EPR line intensity. The intensity should increase with the temperature decrease 
due to the fact that the ground state becomes more populated as the tem- 
perature is lowered. The functional dependence of the EPR line intensity 7, on 
temperature is given in [1]: 


ep | —5.5(5+ 0 |/k,T 
a= 


S 


Pai PRAĆ WOLF (3) 
) (25+ 1)exp | --JS(S+ 1)/k; r| 


Fitting the experimental data on the intensity changes with temperature, 
the value of J = 5.0(5) cm '* was obtained for GdBa,Cu;O> _; [9] and J = 20 
to 30 cm"! for Gdy.sNdoy sBa,Cu30>-; [10]. 

Apart from the EPR spectra corresponding to Cu”" ions one also detects 
the spectra of rare earth ions when they are the Kramers ones. The spectrum 
of Tm** in TmBa,Cu;0,_; was reported in [11] and is presented in 
Fig. 7 (low field portion of the spectrum). The analysis of this signal, carried out 
in the theory of Lea et al. [14] leads to the conclusion that the signal is due 
to the ['„ doublet resulting from 4f'*(7F>,,) ground state split by a cubic 
crystalline field. The measured g — factor could be explained assuming that 
the phenomenon of the reduction of the orbital moments is essential. 

The EPR signal of Nd** in NdBa,Cu;O07_; was reported in [12]. From 
the measured g-factors and the theory of Lea et al. it was deduced that the 
observed signal is due to the T$* component of the ground state *Is,, split 
by a crystal field of a cubic symmetry. This doublet can be built from the 
|+Ł1/2>,|+7/2), and |+9/2) wavefunctions as follows: 


2+ 


a,|£9/2)+a,|+1/2> +a,|+7/2) (4) 
with a, = 0.75, a, = —0.60 and a, = 0.25. The values of these parameters im- 
ply that the type of coordination of Nd”* is that of a four coordinated 
tetrahedral or eight — coordinated cubic one, and that the ground state 


zero-field splitting is dominated by the six degree crystal field terms. 

Yet another interesting case of the observation of an EPR spectrum of 
Tm? * in the tetragonal phase of the (Lu, Tm) — Ba — Cu —O system has recent- 
ly been reported in [13]. The EPR signal of this ion is presented in Fig. 10 
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Fig. 10. The EPR signal for the tetragonal phase of (Lu, Tm)—Ba—Cu—O 


in the low-field part of the spectrum. The observed EPR line could be described 
by the following spin-Hamiltonian: 


= A — 


x =HUpSgB+IAS (5) 


with the effective electron spin S= 1/2 and the nuclear spin I = 1/2. The 
spin-Hamiltonian parameters determined directly from the observed spectrum 
were: g = 3.12(1) and A = —952 MHz. These values are somewhat smaller 
than those found experimentally for Tm** in CaF, and those which could be 
deduced from the crystal field theory of Lea et al. applied to the *F>,, ground 
state of Tm?*. This fact again supports the idea that the reduction of the 
orbital angular momentum, lowering both the value of g and hyperfine struc- 
ture parameters, is essential in the discussed compounds. Physically, such redu- 
ction is due to the delocalization of magnetic electrons into neighbouring or- 
bitals of ligands, or in the other words, is due to rather pronounced effects 
caused by overlap and covalency. 


EPR MEASUREMENTS OF CONDUCTIVITY 


A generalized molecular crystal model (GMCM) developed by Triberis 
and Friedman [17] is used to study the small-polaron hopping motion in 
disordered systems such as discussed here. 

The treatment is based on a generalized „hopping-model” Hamiltonian of 
the following form: 


<m|H|n) = śm|HG+V|n> = E; in OO nan + LM] V |n). (6) 


Here ćm|V |n> describes the overlap interaction, |n> = li, (n,)> are the eigen- 
states of H, and H is the zeroth-order Hamiltonian of Holstein the eigenvalues 
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of which are: 
|| 
EF e;(0) — E,() + ), ho, (n+>). (7) 
k 


where (n,) represents the totality of the vibrational quantum numbers 


(..., n,,...) for the occupation at r,, 


E,(0 = +7 A? 2Mof 8) 


is the small-polaron binding energy, and e;(o) is the local electronic energy. 
N is the number of lattice sites and 4, the electron-lattice interaction 
energy. 

The GMCM has essentially two advantages: 

(i) the local electronic energy e,(o) is site-dependent, 

(ii) the electron-lattice interaction parameter A, and the binding energe 
E,(i) are also site-dependent. 

According to Triberis and Friedman for the high-temperature regime the 
conductivity varies with temperature as: 


o = exp[ —(T,/T)**] (9) 
with 

1, =8.5N] a" "/kgN;, (10) 
where N, is the concentration of sites, « " is the spacial extent of electronic 
wave function localized at a single site, and N, is the density of states. 

In the low temperature range the conductivity varies as: 

o = exp[ —(T„/T)"*] (11) 

where 


T, = 6.49a2/k5N,. (12) 


Direct measurements of the resistance for Pro,.sREg sBaaCu;O0, (where 
RE = Tm, Lu, Er, Yb) shown in Fig. 11 supported unequivocally the applica- 
bility of the discussed model [13]. 

As it is well known the increase of the conductance affects strongly the EPR 
line intensities. This phenomenon is mainly caused by the skin effect which 
reduces an effective sample volume for penetration by microwave radiation. 
The skin depth 5 for microwave penetration is inversely proportional to the 
sample conductivity: 


5 = (2/u,wo)"?. (13) 


5 — Szczecińskie Roczniki Naukowe t. IX, z. I 
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Fig. 11. Temperature dependence of the relative resistance for Pr, „RE, „Ba,Cu,O, 
(RE = Lu, Yb, Tm, Er) 


Semiconducting behaviour was observed for samples with large oxygen deficie- 
ncy. The relative intensity changes of the EPR line for semiconductors can be 
described by a formula given by Godlewski et al. [4]: 


ysg_ 2exp(—w)+[1—exp(—2w)]w"" 
jad [1+exp(—w)]” 


where: w = d/o (here d is the thickness of a sample). This expression was used, 
for example, to study the mechanism of the conduction for GdBa,Cu;O0-_5; in 
its tetragonal phase [5]. Table 1 presents the measured and calculated from 
(13) and (14) parameters for this system. 


(14) 


Table 1. Measured and calculated parameters appearing in (13) 
and (14) 


Temperature Relative % Ó o 
[K] intensity [mm] [Q7' cm ''] 
30 l — — — 


The conductivity a was calculated from the measured intensity ratio I/I, 
and from equation (14). Table 2 presents the ratios lno;/lna,, both experi- 
mental and those implied from the model of Triberis and Friedman, for dif- 
ferent temperatures T; ź T,. 

Good agreement with the theory of Triberis and Friedman is observed. 
This agreement is partially due to the fact that the non-stoichiometry 
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Table 2. Experimental and calculated rations of lno;/lno;, 


Ina5, /IN0yp 1.27(10) 117 


In as, /Ina,,, 1.12(5) 1.12 
CC 1.09 (5) 1.06 


of GdBa,Cu30-_; introduces disorder in the system. Taking the random 
arrangement of the constituents and the above results into account one 
is entitled to claim that the conduction process in this material can be 
attributed to small-polaron hopping motion between copper 1ions in different 
valence states. 

These results seem also to be in good compliance with the idea that 
the presence of small-polarons prevents the appearance of superconduc- 
tivity. 


CONCLUSIONS 


1. The static magnetic measurements give fairly exact information 
about the critical temperature for high 7, compounds. For the compounds 
with large deficiency of oxygen (ie. for the compounds in tetragonal 
phase) the static magnetic measurements revealed the existence of many 
unusual types of magnetic phase transitions and long-range magnetic 
orderings. 

2. EPR studies of Cu?* ions indicate the fact that the local symmetry 
of these ions is orthorhombic for superconducting samples and axial for the 
samples with large oxygen deficiency. 

3. At low temperatures, usually below 25 K, one is able to detect the EPR 
signal resulting from exchange interactions between Cu** pairs. The tempera- 
ture dependence of the signal intensity allows the estimation of the value of 
exchange integrals. 

4. The Kramers rare-earth ions are also frequently detected by EPR even 
at high temperatures. The common feature resulting from such studies is 
the reduction of the angular momentum hyperfine parameters. This phenome- 
non suggests that the rare-earth ions in the studied systems are bound with 
its ligands not purely ionically but there is a contribution from overlap and 
covalency effects. 

5. The EPR studies of the temperature dependence of the line inten- 
sity carried out at high temperatures for the samples with large oxygen 
deficiency can be explained in terms of the skin effect. The deduced values 
of temperature variations of the conductivity are in good agreement with 
the theory of Triberis and Friedman suggesting that the small-polaron 
hopping mechanism is responsible for the conduction process in these 
samples. 
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NIKOS GUSKOS, JERZY KURIATA, VLASIS LIKODIMOS, LUCJAN SADŁOWSKI, MARIAN WABIA, FRANCISZEK 
LEMBICZ, WITOLD DUŁAK 


BADANIA STATYCZNE MAGNETYCZNE 
I BADANIA 
ELEKTRONOWYM REZONANSEM PARAMAGNETYCZNYM 
WYSOKOTEMPERATUROWYCH NADPRZEWODNIKÓW 
I ZWIĄZKÓW POKREWNYCH 


Streszczenie 


Praca stanowi przegląd możliwych zastosowań statycznych pomiarów magnetycznych i,elek- 
tronowego rezonansu paramagnetycznego w badaniu wysokotemperaturowych nadprzewodni ków 
i związków pokrewnych. Podano i przedyskutowano szczególne przypadki tych zastosowań. Wy- 
kazano, że pomiary magnetyczne ujawniają istnienie różnych magnetycznych przejść fazowych, 
natomiast pomiary EPR dają unikalną informację o oddziaływaniu szeregu centrów paramag- 
netycznych z ich krystalicznym otoczeniem. 


HHKOC TYCKOC, E)KbBI KYPATA, BJIEACHC JIMKOJJMMOC, JIBYLIMAH CAJJIOBCKH, MAPHAH BA5B4, 
©OPAHLMLIEK JIBSMBHU, BATOJIBĄ NYJIAK 


CTATHUECKME MATHHTHPBIE MUCCJIEIOBAHHA 
M HCCJENROBAHNH4A 
SJIEKTPOHHbBIM IIAPAMAIHHTHbBIM PE3OHAHCOM 
BBICOKOTEMIIEPATYPHbBIX CBEPXIIPOBONHMHKOB 
HM POHCTBEHHbIX COEHJMHEXMM 


Pe3roMe 


PaOoTa aBJIaeTca 0630p0M BO3MOXHbIX IIDAMeH€HHA MAaTHUTHbIX M DIIP u3MepeHuń npu 
KCCJIEĄOBAHHM CBEPXIIDOBOJHHKOB HU DPOĄCTBEHHbIX COEJJMHCHHi. OOcyXK1aeTCcA OCOÓBIe CJIYMAH 
TAKHX IIDHUMEHEHHA. IIpeHCTABJIEHO, UTO MaTHMTHbIe U3MEPEHHA IIDOABJIAKOT IIDACYTCTBHE MAaT- 
HMTHbIX (hAa3OBbIX IIepeXO/1OB, a u3MepeHHa DIIP NqaroT YHHKAJIbHylo UHQOPMALNIO O B3AHMOJCH- 
CTBHH pANA IIADAMATHHTHbIX IIEHTDOB C HX KPHCTAJIJIMUECKAMH OKPECTHOCTAMA. 
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Summary: The aim of this paper is to develop the theoretical basis of a deci- 
sion choice under risk conditions and to present a batch of the computer program- 
mes with the help of which decision making is easier. The theoretical model will be 
constructed based on the game theory methods. 


Key words: stochastic planning — game theory — linear programming. 


In the present paper we consider a decision problem (%, %, 7) where % is 
a set of the decisions available on the information level of a decision maker, 
R is a set of the constraints and 9 is a plan of the decision maker preferences. 
The decisions x e Ż satisfying the constraints form a subset X of the feasible 
decisions. The plan ? determines the optimum decision from X. One assumes 
that the plan ? is constructed with the aid of one or a few functionals /,, ..., f, 
defined on %. The mapping f =(/f,, ..., f,) transforms a set Ż into the criterion 
space K = f(%) c E*. The plan of preferences can consist of the determining in 
K (or E*) an ordering relation go, which generates in % the relation 


ć no worse than u <+f(ć)ofln), ć,neX 


permitting to compare the decisions with regard to their usability for a decision 
maker. Usually, the appropriate determination of such a decision is difficult 
therefore one makes the additional assumption that each of the decisions xe % 
can be characterized by a system of numerical parameters (x,,..., x„). Then 
the set 4 becomes a subset of a space E” and the functionals f, become the 
functions of n variables. If it is possible to present a set % as a system of 
equations or inequalities determined on % (or E”), then the determination of 
the optimum decision is the question of the mathematical multi-decision pro- 
gramming which is well described theoretically and also numerically (e.g. in 
[1, 2; 12]). The application of the decision models of this form is bounded as 
they do not take into consideration the occurrence of the random factors. 
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Therefore, we will further assume that a decision model has been constructed 
in a form as above but at least some of its parameters are the random variables 
determined on some probability space (2, M, P). Such a model will be called 
the stochastic decision model. In particular it may be a plan of the enterprise 
activity. Constructing of such plans is especially important for a farm enterprise 
because in this branch of production the random factors play especially impor- 
tant part. Further considerations will, first of all, be concerned with such type 
of enterprises. 

Let C = (Cz, C,) be a random vector and its components be all the random 
variables of a model and let A = (Ap, 4/,) be a set of the values of a variable C. 
The vector Cą contains the components of a vector C occurring in the const- 
raints and the vector C, contains the components occurring in the functionals 
J,, Analogously Ap, A, contain all the possible realizations of the vectors 
Cz, C,, respectively. Replacing the vector C in a model by its realization Ce A 
we obtain a model determined by a decision set X, and numerical functionals 
Ji. Fixing a plan 9 we can find an optimum decision x„e X,. However, under 
other realizations it can be infeasible or very disadvantageous. This first dif- 
ficulty can be avoided by taking as a decision set the product 

Ze Fw, 
ceAR 
under condition that it is non-empty. To reduce the second difficulty we can 
use the additional information concerning the either possible realizations or the 
distributions of a variable C. 

In agriculture one employs various methods of a stochastic planning for 
about 50 years. The papers [5, 7] contain a comprehensive and critical review 
of the bibliography concerning this subject. The aim of this paper is to develop 
the theoretical basis of a decision choice under risk conditions and to present 
a batch of the computer programmes, by the help of which a decision making is 
easier. The theoretical model will be constructed basing on the game theory 
methods. The necessary definitions and theorems of this theory can be found in 
the papers [3, 11]. There are a number of papers describing the applications of 
game theory in agriculture. In Poland these investigations were made by a few 
authors only. The best results which could be trained into a planning practice 
have been obtained by A. Miklewski. His papers [10, 9] contain also a comp- 
rehensive review of literature. The using of the recent methods especially 
multi-decision and nonlinear programming allow us to obtain the more general 
results from a theoretical as well as numerical point of view. 

In the farm enterprises the radomness of the model parameters arises as 
a result of the unpredictable alterations of the natural and the market con- 
ditions. The set of factors creating the distributions of the model random 
variables will be called a state of nature. The subsets of a set A or the families of 
distributions will compose the nature strategies sets. The subsets of the set 
% will form the planist strategies sets. With the aid of the functionals f, 
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the profit function and the loss function are constructed in a different way. 
They are either the random variables depending on a variable C, immediately 
or the statistics depending on tis realization (samples). They have to express the 
basic decision maker preferences but they cannot be too complicated because 
of the difficulties in a collection of the necessary observation materials, as well 
as of the game solution. These general rules will be further defined according to 
the undertaken assumptions about the information being in receipt of a deci- 
sion maker and concerning either the possible realizations or the distributions 
of the variable C. 

At first we assume that a planner does not have any information about the 
distribution of the variable C at his disposal but on the base of the statistic data 
or of his own observations he can determine certain subset A” = (Ak, 4) 
of a set A. If a decision maker does not take into consideration any risk to 
exceed the constraints then he can replace a random variable Cz by its expec- 
tation E(C3) which based either on data from a set A or on the additional 
non-model information is prognosed for a time of a plan realization. The set A” 
being finite then E(C4) = c is a mean realization. In this way a decision set 
X, independent of hazard is obtained. If such a solution is considered as being 
too risky the set X is used. One calculates a random profit function W(X, C,) 
for each of the realizations c,e Ay. Now one can consider a game (X(or 
X), Ay, W). If the functionals f, are linear then as W (or L) two-linear func- 
tion x"Oe, can be taken. In this case the use of other loss function 
(for example a square one) would complicate a question too much without 
any special profits. To construct this function in the multi decision models 
it is most convenient to apply a method either of the limenial constraints 
or of the utility function. If the sets X and A; are compact and a profit 
function is continuous on X x A, then there exists an upper and a lower value 
of the game v, and v, > v,. The most safe planner strategy is max-min strategy 
x, assuring him under each nature strategy a payment v, at least. Namely 
we have 


max min W(x, c,) = W(xy, ej) = v,. (MAX) 


xEX cyqEAf 
Min-max nature strategy c, determined by 


min max W(x, c,) = W(x, cy) = vj. (MIN) 
cjeAf xE€X 
makes impossible that a planner realizes the profit superior than v, i.e, 
he minimizes his loss. The loss function L can be treated as a nature profit. 
It is its "profit" and a loss for a planner. A planner, operating rationally, 
will try to limit a nature profit at minimum using for this aim a min-max 
strategy realizing 


min max L(x, c,) = L(x, cy) = 04. (S) 
xEX cyq€EAf 
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If the sets X and Ay are compact and convex, for example if they are the 
polyhedrons (for a finite Ay it may be conv Ay) and a salary function 
is two-linear then v;, = v, = v and the max-min strategy is an optimum planner 
strategy. 

In general, there is not any possibility to determine all nature strategies, 
and even if it would be possible and if an optimum strategy would exist 
then the necessity of taking into consideration an extremaly unfavorable 
but few probable nature strategies would cause a payment of very small value 
that would be difficult to accept by a decision maker. The first possibility 
can easily be explained by the fact that nature is not a rational player but 
it uses some fixed (but unknown) random actions to choose its states, and 
thereby also the strategies of a game determined by a planner. Thus, there 
is no reason to assume that in every such a game he will apply his 
optimum strategy. The second possibility avails itself of the fact that 
a planner has almost always some information which of the possible 
nature strategies are so little probable that they can be abandoned with 
a small risk. 

Now we shall more consider in greater detail the case when 4? is a finite 
set, 1.e. 


C' = (Ci, C3,..., C,) 


C (OC: „w Ch 
If all constraints are linear, ie. if they are of the form 


a,' 


ł 


xsb, (G=1,2,..., PD) X. 0, (*) 
then the set X is a convex polyhedron determined by the inequalities 
PEM =(= 124.10); (ze) 


arisen from (*) when we replace the random coefficients a, and b, by 
their realization determined by ch e Ap (we assume that the random constraints 
are of inequality form). If limenial constraints method is used, then to 
the inequalities (+*) one should add the constraints f;/(x) > (<)y,(i=1,..., k) 
suitable for all strategies cy. The value of a distinguished functional i.e. 
W(x, e,) = r;(x) = fJ(x) is a salary function. Thus, for i 4 s the functionals 
f/ influence only indirectly a value of a payment by the limiting of a decision 
set. In the utility function method a linear combination r;(x) =) a,f; (x) 
with the non-negative coefficients a, being a two-linear function of the 
vectors x and c,, is the salary function. Under above assumptions a game 
(X, Ay, W) =<(X, cj r, has in general no solutions and thereby the 
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optimum strategy of a planner does not exist. To find a max-min strategy one 
should in accordance to (MAX) 

(a) to maximize v>0 
under the constraints: 

(b) r;(xjżv for j=1,2,...,r; x200, 

(c) Bx<$P, x->0; 
The constraints («*) are belonging to the constraints (c) determining the 
set X and also the limenial constraints. It is a usual question of linear 
programming. 

Instead of a salary function r,(x) we can use a loss function 


z;(x) = r;—r;(x), 
where 
r, = maxr,(x) =r;(x). 
xeX 
The differences z,(x) are called the regrets (of Savage) because they deter- 
mine the wasted possibilities of a planner who has applied a strategy x instead 


of an optimum strategy x. The minimization of a loss function leads according 
to (S) to the searching of 


Bow 
and of a strategy x, realizing this min-max. It leads again to the question of 
linear programming 

(a,) minimize u>0 
under the constraints: 

(Gz kież>5 fot Jedszanąsky , K20: 
and under the remained technical constraints: 

(c,) 8x <f, x>0, 
contained also the eventual limenial constraints. Obtained as a result of a sol- 
ving of this problem, the planner strategy raises in general a salary for the mass 
of nature strategies, however it admits also for some of them a salary lower 
than v,. 

The strategy can be reduced either on the basis of a non-model informa- 
tions about the chances of the particular realizations or on their estimation in 
view of the data contained in a model. One of the possibilities is the ordering of 
the realizations with regard to a rule fixed by a decision maker, for example, 
with regard to a "distance" from the realization of a mean in a metric sense or 
with regard to the value of a distinguished functional. The rejection of the 
extreme realizations and an application of a max-min strategy to the reduced 
sets gives a number of the solutions, and hence a decision maker can, regarding 
the additional conditions, choose an optimum one. That metod for a small 
number of the observations can be little effective. 
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Now let us assume that a family P,„(ve0) of possible distributions is kaown 
and that an unknown distribution of a variable C belongs to it. In general, the 
information being in receipt of a decision maker allows to state or to estimate 
only some parameters of these distributions. If only their expectations E (C) 
are knowi, taking A = (E,(C): ve©); the problem can be reduced to the pre- 
vious one. It does not need to use either analytic or numeric methods. If the 
covariant matrix S$, (ve ©) of a variable C is also known then a decision maker 
has much broader possibilities. Often, the variables in both vectors Cę and C 
are weakly reciprocally correlated, which permits him to treat them separately. 
If it may be assumed that the variables C; have a normal distribution then 
a decision set can be constructed with the aid of the constraints made on the 
probabilities of an inequality set realization (a Charnes — Coopers method). 
However, it needs the non-linear constraints and causes that the problem beco- 
mes much more complicated especialy for the greater number of the random 
constraints. Considering the present possibilities of the collection, storing and 
transformation of the information we think that the using only the expectations 
E,(CĄR) or the information about the realizations Aą in the way mentioned 
above is sufficient to construct a decision set. The information about the cova- 
riants of the variables C, can be used for the construction of square loss 
function. For one random linear functional with the coefficients C, this function 
is of the form 


L(x, v) = ((C,-E,(C,): x)”. 


It is a random variable thus for the game determination we introduce 
a number risk function 


K(x,v5)=E [L(x,v)]=x'S,x, xeX, ve0 


and we consider a game (X, ©, K).S, is a covariant matrix of a vector C g for 
the distribution v and x" S,x is a variance of a random functional f = )C;x;,, 
K is treated as a loss function. Usually it is possible to estimate the covariant 
matrices for only one distinguished random functional. Then one can use a li- 
menial constraints method. It consists now in addition to the constraints of 
a set X, the limenial constraints of the form E,(f, > (<)y,i=1,2,...,k,and 
they together determine a decision set X, as well as a game 4X, ©,K). To 
determine a mini-max strategy x, the calculation 
min max x" S,„x = xySXy (=) 
Xp 0 
The planning practice realized up to now has showed that usually we lack 
either the observation materials or the statistic materials to estimate the cova- 
riant matrices for different hypothaetic distributions v and we content ourselves 
with only one such a matrix S$. Thus the formula (*) passes into 
minx"Sx = XxySXy, XMEX5. (se) 
Xp 
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The question is reduced now to the minimization of a variance of a distin- 
guished functional f, on a set X,. If for f,, i * s the constraints are fixed and 
a coefficient y, is alternate then the formula (**) determines a well known 
method VE of stochastic programming. 

If, in the enterprise, there exists a possibility to construct a mathematical 
model of the form such as above then the batch of computer programs based 
on a described method gives the road for possibilities of a decision choice 
under risk. It makes possible: 

I. The usage of linear methods i.e.: 

1. A determination of a set X, or X (according to the decision maker 
wishes) under the realizations Az either observed in the enterprise or obtained 
on the basis of non-models data (*) as well as an application of the mul- 
ti-decisions optimization methods, for example a limenial constraints method, 
an utility function method, a lexicographic method, a STEM method. 

2. A determination of a max-min strategy and of the Savage's strategy (for 
a set X) using additionally either a limenial constraints method or an utility 
function with an eventual limiting of the strategies. 

II. A minimization of a risk with the aid of a risk function K. A set X, 
remains without change. If a covariant matrix for a distinguished functional f is 
known, then the programme uses a VE method to review non-dominated 
solutions for a couple (E,, o,) with a variant choice of limenial constraints for 
other functionals. Based on our own programmes — a linear and that of 
a variance minimization with an application of a modificated algorithm 
of Frank — Wolfe [4, 8], we have constructed our batch. The programmes are 
written in Pascal. 

The possibilities of the batch will be presented on the following simple 
example concerning the vegetables production plan of a farm enterprise of 
a middle size and the optimum one with regard to the climatic and the soil 
conditions occurring in the valley territories of Poland (a modified model 


from [6]). 
The cultivated vegetables and their acreages: wheat — x;,, rye — x,, 
rape — x, potates — x,4, corn — x;. 


A decision maker is going to optimize a size of the production taking into 
consideration the random changes of a capacity and a non-random deficit of an 
organic soil substance. 


Constraints 


l Xył Xą+ X4+  X4+  Xs=96 balance of grounds 

2 X X <60 _ limitation of grain areas 
3. Xą < 14  [imitation of rape areas 
4.  8.2x,+8.0x, +8. 2x;+32. 1x,+8.1x; £ 1640 balance of a work 

5 0.8x; + 4.8x, ++ 2.0x; £ 330 work top I 

6.  4.5x, +4.7x+1.9x3 x290 work top II 

7.  D,x,+D,x,+D3x3 > 2170 random balance of straw 
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Objective function — a total production 
Y = Ć,% + C,x, +C2Xx3 +C,x,+C;sxs w TZICO) 
Deficit of an organic soil substance 
R= —1.5x, — 1.5x, — 1.5x3 —3x4—3Xx5 = Ja(X), 
C= (SE Cz, Cz, C, C;, D4, D,, D3), 
C,=(C,, Cz, Cz, (4, €5),  Cęz=(D4, D;, D3). 
Observed realizations of a vector C 


Cą Ca Cz (of Cz d, d, dą 


C;: Cz 29 63 35 
47 172 31 
54 60 39 
Means 54 64 34 
Covariant matrix of a vector C, 
96.04 52.65 0 0 0 
52.65 46.24 50.38 0 0 
S: 0 50.38 42.25 0 0 
0 0 0 150.68 0 
0 0 0 0 134.65 


In a computer memory the coefficients of constraints 1 —6 and of a func- 
tional R, as well as the vectors matrices C,, Cą (realizations) and a covariant 
matrix S$ of a vector C, are well-remembered. From it one can easily observe 
that the 5-th planner strategy giving the best production result (against the 5-th 
nature strategy) is unfeasible against the 4-th nature strategy because the 
left side of the 4-th constraint equal to 1847.3 is smaller than 2170, and on 
the other hand, a very profitable 3-th planner strategy (against the 3-th 
nature strategy) giving the production Y =4031.05 and R= —197.92 
gives a very low production value (Y = 3576.42) against the 4-th nature 
strategy. 

The extremum values Y and R on the set X, and X do not differ what 
shows Table 2. One can recognize that the values of R are to great, however 
changing a coefficient, of the utility function a decision maker can obtain more 
advantageouns (according to his preferences) relation between Y and R 
(a programme gives such a possibility). 
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Table 1. The sets X, suited to the nature strategies (C,) as well as the planner strategies realizing 
a maximum of the functionals f, (without any constraints for a functional R) 


ea |a| | | pt rem 


| | = 
ż | < 
3 0 = 
4 8.2 < 
|-> 0 < 
6 4.5 < 
| 53 60 34 0 0 2 
2 57 65 31 0 0 2 
| 3 Sy 63 35 0 0 2 
4 47 72 31 0 0 2 
3 54 60 39 0 0 2 


Objective function | 


1 30 38 
2 5 32 h 39 
A 40 31 35 
4 32 36 40 
5 37 29,5 " 45.5 
| The optimum planner strategy | 


3825.04 
3829.50 
4031.05 
3959.41 
4189.37 


3861.69 
3846.6 
3324.0 
3324.0 
3702.57 
3702.57 
3550.67 
3509.87 


(e) 
5 


ZJ A A WJ II I R 


Objective 


function Y: 
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Table 3. A max-min strategy and a strategy of Savage for a set X under a. limenial constraint 
R > —200 


Table 4. The values of functionals Y and R for max-min and Savage strategies and all nature 
strategies 


Planner Nature strategies 
strategies - m="L 
| 


3138.29 3757.83 3931.32 3738.29 3863.28 
3770.57 3795.61 3979.65 3716.38 3904.02 


4 l 
Table 5. A model max-min for a utility function u=—: e: "R=0.8:Y+0.2:R 


| 


5 


cą 
1 


[50 
OCOOCOANDOOCH 
O — 
OOOpWCO- 
C© — 
©.* Oe.v.v ve 
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Objective 
function 
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Table 6. A model of Savage 


- 
14 
1640 
330 
290 
2170 
2170 
2170 
3009.69 
3023.00 
3185.33 
3103.89 
3262.51 


|so) 
OOOCOANOCOOCO- 
Go — 
OOONE%WGOOHL 
CO 


sa) 

* O 

= 
© 


Table 7. The optimum planner strategies for particular nature strategies 


3009.69 
3023.00 


3185.33 
3103.89 
3262.51 
3060.68 


Table 8. The max-min and Savage strategies for a model from tables 5 and 6, respec- 
tively 


[==] | LE [ETL [= 


6 — Szczecińskie Roczniki Naukowe t. IX, z. 1 
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Table 9. The values of functionals Y, R for the max-min and Savage strategies and for all nature 
strategies 


Nature, strategies 
Planner pa z aa 


3194.20 3799.17 3935.24 3794.20 3993.20 214.99 
3176.69 ZB 28 3896.37 3801.57 3992.85 4 


In part 2, for an arbitrarly chosen constraint for R (for example R > —200), 
computer makes a graph, presenting relation between standard deviation of 
a distinguished functional f, and a minimum value of an expectation of a varia- 
ble Y (ie. under a constraint that Y ża). These graphs for R> —180, 
R > —200, R> —220 are shown in Fig. 1. 

The ought acreages for a given R under the most profitable for a decision 
maker, relation between E(f,) and o, can be read from the graphs in 
Figs. 2, 3 and 4. 


ST. DEU 


38 38 37 38 389 
<X 206) 


PROOUCT ION 


Fig. 1. Relation between standard deviation of a distinguished functional f, and a minimum value of 
an expectation of a variable Y, for R> —180, R> —200, R> —220 
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Fig. 2. Acreages for R = — 180 under the most profitable for a decision maker relation between E( f,) and o, 
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Fig. 3. Acreages for R = — 200 under the most profitable for a decision maker relation between E( f;,) and o, 
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R=-220 
U LJ U 
46 |-_. 
whaat 
p potatos 
| reaterar iti taRtAŁG, "| 
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rape e.) 
16 i . a A TOWNEC PI PARKA śezwoini te NO dezti ANA me 
* i eorn 
k I 
+ i 
Led 
x. 
h A 
4, 
"A - 
o % ? 
|) PISZ AAA gda Pw AA OZ AFC EOT PJ SA O waz AA ' ' ł U U =" manietyniaw yte .- R+++++"TYUE"""""" res 
= Z i | i 1 1 [I l 
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PROOUCTION 


Fig. 4. Acreages for R = —220 under the most profitable for a decision maker relation between 
E(/, and o, 


In view of a linear part Fig. 1 a decision maker is fixing a feasable range of 
a variability of the parameters E(Y) and R and basing on a second part he is 
choosing a feasible solution minimizing a risk. 


JÓZEF KOPEĆ, ZBIGNIEW GUZIENIUK 


STOCHASTYCZNE PLANOWANIE DZIAŁALNOŚCI 
PRZEDSIĘBIORSTWA 


Streszczenie 


Celem pracy jest rozszerzenie podstawy teoretycznej wyboru decyzji w warunkach ryzyka oraz 
przedstawienie pakietu programów komputerowych mogących ułatwić jej podejmowanie. Model 
teoretyczny został skonstruowany na podstawie metod teorii gier. 
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10330 KOII3LIb, 356HMTHEB TY3EHIOK 


CTOXACTHMUECKOE IJIAHMPOBAHME J[IEATEJJBGHOCTM 
HIPETNPH ATMIA 


Pe3roMe 


LleJibro pa6OTbl ABJI4eTCA paCIiMpeHue TeOpETHUUECKHX OCHOB IIDHHATH4A pELLIEHHIA B YCJIOBHAX 
pucka. IIpeqcTaB.IeHbl KOMIIBFKOTEDHbIe IIDOTDAMMBI, OÓJIETUAIOINHE IDHHATME TAKOTO pPELIIEHHA. 
TeopeTuueckaa MOJIEJIb IOCTPOEHA Ha OCHOBEe METOĄOB TeOPHM HID. 
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TRACKING SYSTEMS IN H, SPACES 


Institute of Control Engineering, Technical University of Szczecin, 
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Summary. The non-linear tracking systems are considered in the paper. The 
idea of optimality and strong optimality of tracking systems is given. It is assumed 
that the class of input signals belongs to the space H,. Method of describing function 
is briefly presented. Two theorems for finding the optimal (in the meaning given 
above) non-linear controllers are given. 


Key words: feedback control system — nonlinear — optimization — desc- 
ribing function — frequency methods. 


INTRODUCTION 


Let us consider a feedback control system described by the 
equations 


= G(x(0)) 1) 
x(t) = z(t)—F(y(t)), 


where z(t) denotes the reference input signal and x(t) — the error signal. Opera- 
tions G and F describe the controller and plant respectively. Let operation 
H,„23P: L. > L? depend on plant F and controller G. We take into considera- 
tion only those signals z(t)e L* for which w(tjel?* and z(t) = P(w(t)) exist. 
Let Z(P, m) be a set 


Z(P, m)=fz: z,wel”, z =P(w) and |w| <m< o). (2) 


If constant m < o is not determined precisely then instead of Z(P, m) we can 
write Z(P). For simplicity we can assume that P(w) = F(w) for welL?. The 
considered tracking system is shown in Fig. 1. 

Let C be a collection of controllers which can be used in our control 
system. 
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Fig. 1. The considered tracking system 


By xg(t) let us denote an error signal in the control system (1) with a cont- 
roller GeC, obtained for some z(t)e Z(F). 

Let us suppose that for z(t) e Z(F) there exists a solution xq(t)e Z(F) of the 
system of equations (1). 


Definition 1 


A controller GoeC is said to be optimal for input signals z(t) e Z(F) (op- 
timal for a class Z(F)) if for every GeC there is 


sup ||xG.l|| — sup |xGll £0. (3) 
zeZ(F) zeZ(F) 


The control system (1) with the controller Gy is said to be optimal for the class 
of input signals Z(F)(Z(F, m)). 
Definition 2 

A controller GyeC is said to be strongly optimal for input signals 


z(t)eZ(F) (strongly optimal for a class Z(F)), if for every GeC there is 


sup (lxGll — lżcgll) < 0. (4) 
zeZ(F) 
The control system (1) with the controller G, is said to be strongly optimal for 
the class of input signals Z(F) (Z(F, m)). 
Formulas (3) and (4) can be rewritten in the form 


iE= min (sup 4 xl: weH,, ||w|| £ m < 0o)), (5) 
I, = supfmin |x|: weH,, ||w|| < m < ooj. (6) 
GeC 


The equivalent standard problem in H _-control theory is obtained by defining 


(see [47) M _ k =H H 


y = G(x). 
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Pe") 
The standard block diagram is shown in Fig. 2 where G= Ę A 


Fig. 2. The standard block diagram 


The following example illustrates the idea of the optimality and strong 
optimality of feedback control systems. Let m in the definition of Z(F, m) 
be precisely determined. Using definition 1 we can answer the question: 

What is the least portion of energy that should be delivered to the system in 
order to track any reference signal z(t)eZ(F, m)? 

If we are interested in saving the energy, we should employ definition 2. 

Let us assume that operation F is a linear one given by the formula 


F(y(0)) = | y(t—1)df (2), 


where f (1) is a bounded variation function or an operation given by the trans- 
fer function F(s). On nonlinear operation G we impose additionally the con- 
dition y(t) = O for x(t) = 0. According to the above assumptions the equations 
(1) can be written in the form 


y(t) = G(x(t)) 


i (7) 
x(t) = z(t) — | y(t—c)df (2). 
0 


DESCRIBING FUNCTION METHOD 


Now we present briefly the ideal filter method described for instance 
in [5]. Let us take into account the equation (7) and assume that for z(t)eL* 
there exists exactly one solution to equation (7), which belongs to IL. 
Let z(t) = z,(t)+z,(t) be a sum of almost periodic function (belonging to 
the Banach space of a.p.f.) and a function from the space L? such that 
z,(t) = —z,(t) for t < 0. A solution of the equation (7) takes the similar form 
x(t) = x„(t)+x,(t), where x,(t)ea.p.f. and x„(tjel? with x,(t) = — x,(t) for 
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t<0. The equation (7) can be expressed as a system of the following two 
equations 


x,(t) = z,(t)— j G(x,(t-n)dfo), t>0, (8) 
a 


i 


x(t) = —|[G(x,+ x.) — G(x,)](t—1)d/ (2) + z, (t) 


0 
+[IGOJIC-IYO0, t>0. 0) 


It can be proved that for z, (t)e L* and z,(t)ea.p.f. all solution of the equation 
(9) belongs to I. 

Now we consider the equation (8) in the space of almost periodic (or perio- 
dic) functions. It can be rewritten in the form 


[e ©) 


x(t) = z(t)— | G(x(t—1))df(r), te(—0, 0). (10) 
0 


Let us suppose that z(t) is a periodic signal with period T = 2n/o. We look 
for solutions of equation (10) with the same period. Let us express functions 
z(t), x(t), y(t) as Fourier series 


zb="Nz;e" = „AE eĘ zM=RIe". (11) 


Putting extensions (11) into the equation (10) we obtain 
X, x FOK O Ra. PLZ, ie; (12) 
where 
jp: , JURA | ! 
JalXos X: X...) == |ylhe *"*dt => | [GQQ,xe**)|ltje "dt. (13) 
R 0 T; k 
The system of equations (12) regarded as an operator equation in the space 
of infinite 1” sequences. Introducing 
PZŻ Z Z. JA AZAXY BIE 
V=UoVvV2---), F'y' =4F()yo, F(jlo)y,, F(żjo)y...), 
G'x' m tolxo> X1; X2 z; J1(X0> Xi Xz-.-), JalXo. Xi X2 *> ). „s 
we can rewrite equation (12) in the form 
x =z —F'G'x'. (14) 


Let us denote by P a projector selecting from the system of equations 
(12) the equations with subscript 1 and by (I—P) a projector selecting 
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from equations (12) all equations except for the equations with subscript 1. 
The equations (14) can be rewritten as a system of following two equations 
Px' = Pz' —PF'G'x', (15) 

(I—P)x' =(I-P)z' —(I—P)F'G'x'. (16) 


It can be proved that if the equation (15) possesses a unique solution, then 
so does the equation (16). For example, if operation G satisfies the condition 
with Lipschitz constant « such that 


q =a|F(jo)| < 1, (17) 


then also the following condition holds a|F(jnw)| < 1. A solution of the equa- 
tion (16) takes the form 


(I—P)x' =(I-P)z'—R(Px'+(1—P)z!), (18) 
where R is the resolvent operation. For R the following estimate holds 
|RPx' +(I-P)z))| < — |Px' +(I-P)z!||, (19) 
pot 
where q is given by (17). Substituting solution of the equation (18) into the 
equation (15) we obtain 
Px' = Pz'—PF'G'[Px' + (I — P)z' —R(Px'+(I—P)z')]. (20) 


Assuming that q in the estimate (19) is much less than unity we can take 
into consideration instead of a solution to the equation (20), a solution of the 
approximate equation 


Px! = Pzi— PF'G'[Px'+(I— P)z!] (21) 
OT 
X1 = z, —F(jo)y; (Zo, X1> Z; Z3 ZA (22) 


Assuming that F(jw) 4 0 we can invert operator PF'. Then the equation 
(20) takes the form 


(PF')""Px'—(PF')""Pz'+ PG' [Px'+(I—P)z'] = 
= PG'[Px' +(I—P)z'] —PG'[Px'+(I—P)z' —R(Px'—(I-P)z')|. (23) 
Denoting by 6(x,) the left side of (23) we get 


X1 
F(jo) 


Z1 


$(x,) = a F(j0): 


FIdZy; 813 Za +4.) (24) 


where 


Je ! ! m. 
J+(Z9; X1, Z2 -..) = F | Glzo+x,e/7+x,877*+ > 2;e7% óde "dr. 
0 


n>l 
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For fixed period T = 2n2/w and vector z! = fzy, zą, zą...) the complex 
function y;(z,, x;, Z...) depends on single variable x, and is called the desc- 
ribing function. Equations (21) and (22) following the method of describing 
function become equivalent to the equation 

$(x,)=0. (25) 

In connection with operation G mapping the space L? into itself the notion 
of the describing function can be precised as follows: Let y(t) be a response of 
the system described by operation G to the signal x(t) = N'sinot. The quo- 
tient of the symbolic value of the first harmonic of the output signal y(t) by the 
amplitude of the input is called the describing function and is denoted by 


ż | y(Osin(ot)dr+> [y0) cos (wt)dt 


G(jo, N) = — (26) 


N 


APPLICATION OF THE DESCRIBING FUNCTION METHOD 
TO OPTIMIZATION OF SYSTEMS 


The next two theorems can be proved:. 


Theorem 1 

Let operations F and Ge C map the Banach space I: into itself and, more- 
over, operations G have a uniformly continuous, bounded derivative. If for 
a controller G,eC the expression 


1  óre(ó(jo,N) . —1 _ óre(G(jo, N)) 
re——+——————  im=—+——— 
F . . F . 01 . 
f,(G) = inf |det (jo) AREGCJ (jo) | im(x(jo) 27 
weR l óim(G(jo, N)) 1 5im(G(jo, N)) 


zTYZ5 YA óre(x(jw)) "F(o) 5im(x(jv)) 


attains a non-zero maximum and if for w from the expression (27) the number 


x = sup |F(jno)| 
n=l 
is sufficiently small, then the control system described by equations (7) with the 
controller Go is optimal for the class of signals Z(F, m). 


Proof 


If assumptions of Theorem 1 holds, then for every signal z(t)e Z(F, m) there 
exists a solution x(t) to the system of equations (7) as well as a solution to the 
approximate system of equations. 


t. = G(x(t)) 


t (28) 
x(t) = z(t) — | y(t—c)d(z). 
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Both those solutions belong to Z(F). Let us suppose that the expresion (27) 
attains maximum for some GyeC and wy ER, 0 < N;y. Hence it follows that for 
every GeC and arbitrary oeR, 0 < N these is 


l z, <( l Ę 
+ G,(joy, ——;+ G(jo, w) ! 
(Gaj Śslew. No) |< |(rqoj* U 
Also the inequality holds 
1 -1 ( 1 ję:1 
u —+ Go(jo, N <su Fij* śe. w) 29 
= Ę a)! ou ) KMIWIGNEKCM 5 
N N 
Multiplying both sides of (29) by m > O we obtain 
l z” (7, l =! 
su + Go,(jo, N m | <su — +G(jo, w) m) 30) 
ż: (rs AR ) a ( Fo | | 
N 


Since the Jacobian of operation 1/F+G is non-zero by assumption, then for 
constants 0 < N, < N, either the relation 


1 ża] =, 1 A zo 
(+5) Na] < (+6) wa 


| 1 » =F 1 p zaf | 
(3+5) "N+x > (+6) Ni % 


holds. Because G(0) = 0 then we may assume that the relation (31) holds. The 
inequality (31) implies that for a sequence of numbers 0 << N,<N,<...<N, 
the sequence of functionals 


1 1 
rz 2)? Gdów, N Y| (Gia; „jt Śl. 8) 
(33) 


satisfies the inequalities 


| 1 
(Gia „ŻA v)|< (Gia jaj” Św N))|<- 


for every weR. The functionals (33) describe a family of curves on complex 
plane, which is shown in Fig. 3. 
Let Ź(F, m) denote the set of all periodic signals which belongs to Z(F, m). 
Now let us assume that the expression 


1 zal. 
(FG 2)? Goliw. w) 


(31) 


OT 


(32) 


1 
(ro zj” Uw N ») 


1 
F (65 pt S0e: x.) 
(34) 
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Fig. 3. Family of curves 


on the left side of (30) attains supremum at the point »;, N, for a certain signal 
z,(tjeŹ(F „m) with period 2n/w, (e.g. z,(t) = N,sino,t). Among periodic 
signals belonging to Ź(F, m) one can find a signal with period 21n/0> 
(e.g. z,(t) = N,sinw,t) such that the expresion [(z( jo))/(F (ja))| attains supre- 
mum with respect to w. Let N, and N, be constants such that 


N1z,(J04) N.z,(J02) 
Fa.) Fa) 4 
The inequality N, ż N, holds. Hence, by (34), it follows that 
(rr 1 w | 1 re 
-+Go(jo,, > + G,(jo,, > 
F(jo;) oj 1 N) (FG 0) oj 2 s) 
1 -1 
> 
= j* Śo(i02, NJ) . (36) 
Relations (35) and (36) imply 
N+z,(J04) ( l Zd e 
=" "R ż +G W;, N > 
Fa) |NFGa)" 90% % 
Nąz,(J02) Ę l ! jn 
AZ — + G,(0>, N) 37 
Flo) |NFGwz) ” 9%2 2 SĄ 


On the basis of inequality (37) the expression on the left side of inequality (30) 
can be replaced by the expression 
z. (38) 


1 ARJ A 
(rycz: Śoio. W) Fa) 


sup sup 


zeZ(F,m) oeR 
N 
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Repeating the same argument for the right side of (30) we obtain 


1 282, | z(ja) 
——+0Q N < 
EE ( FE ŚW ) ad i 
< sup sup Gia 5 +G(jo, w) Ea ) (39) 


zeZ(F,m) oeR 
N 


Applying to both sides of (39) the formula |x|| = |[(F"* + G)(F"*"(2))|| we get 
sup |x6.ll< sup [xgl|. (40) 


zeZ(F,m) zeZ(F,m) 


Now let us take into consideration the system of equation (7). If z(t) 
EŹ(F, m), then a solution of the system of equation (7) is also a periodic signal. 
For every N and o there exists a signal ź(t) e Ź(F, m) such that the expansion 
in Fourier series of the output signal from G consists only of the first harmonic 
component. For signals ź(t) the equality holds 


GG. (41) 


Hence it follows that if inequality (40) holds, then also following inequality 
holds 


sup |xeoll< sup |xGll. (42) 


zeZ(F,m) zeZ(F,m) 


Let z(t)e Z(F, m) be a signal belonging to the set of almost periodic func- 
tions. Then the signal F ""(z(t)) in space I? takes the form 


Ą Zą 
R ASS Wh 43 
(6 (jo) F(jo) | ża 


and operation 1/F+G — the form 


RKBE a Żart z JP R a 00: A) (44) 


faz (jw) F (260 ) 


where z,, x, and y, are Fourier series coefficients of signals z, x and y, respec- 
tively. Let the norm of signal (43) equal m, 1.e. 


e|- Glróa|) 77 


Correspondingly, the norm of signal (44) equals 


[e+o)e|-( 


Zn 


F(jno) 


Fo PYn( > X2 X3; ++.) 


> 


Xn 
F(jno) 
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Let us take into account the periodic signal with period 2n/0 
AZ; 
on EyNEK 0, B; OREG (Gu 
5 (Jo) | 


AZ 
F(jo)| 


where a is such that 


ko) - 


1.e. 


j 
- |e|- 


For the periodic signal of the form 


= Ha OO, +.) 
the equality holds 


i ax. | 
—+G |(x)| =|-— ,0,0,..)|=||=+G 
(+ Je" BIE 76 | [(F Je» 
since operation (1/F+G) is one-to-one. For almost every periodic signal 


z'e Ź(F, m), whose frequency equals the frequency of the first harmonic com- 
ponent, such that the formula holds 


[F*+G7"'(F"'2)| = [F*+G7"(F""2))|. 


Hence, it follows that the formula (42) is also justified for almost periodic 
signals. We can write 


sup |xgoll < sup |xcll. (45) 


zeZ(F,m) zeZ(F ,m) 
Theorem 2 


Let operations F and GeC map the space I? into itself and, moreover, 
operations G have a uniformly continuous, bounded derivation. If for a cont- 
roller GeC the expression 


„1_„śre(ódowN) , -1_, óre(ódo,, N) 
z ! F(jwo,)  óre(x(jo,)) F(jo;)  5im(x(jo,)) 
I,(G) = inf |det NAC. M-21 
BIER az l óim(G(joą, N)) a l „Sim(G(jaz, N)) 
A | F(jo;) óre(x(jo,)) F(jvo,)  0im(x(jo,)) 


(46) 


attains non-zero maximum and if for w from the expression (28) the number 
x =sup|F(jno)| 


nźl 
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is sufficiently small, then the control system described by equations (7) with the 
controller G, is strongly optimal for the class of signals Z(F, m). 


Proof 


If assumption of Theorem 2 holds, then for every signal z(t) e Z(F, m) there 
exists a solution x(t) of the system of equations (7) as well as a solution of the 
approximate system of equations (28). Let Ź(F, m) denote, as formerly, the set 
of all periodic signals belonging to Z(F, m) and let zeŹ(F, m). A solution of 
the system of equations (28) is also a periodic signal. Let us suppose that the 
expression (46) attains infimum for some G,. Then we can write 


l zt 1 -1 
sup + G,(jo w) < su Ur + G(jo,, w) 

01,02ER = 1) : * RR F( jo 0) ś 

N N 
Or 
l se l "e 
sup ||->—+)9(0 ) < sup = + 5,(Jo ) „ (47) 
W1,02ER = i ; 01,02ER F( jw; ) : Ź 
N N 


where )7(w) and f,(0) are first coefficients of Fourier series expansions of the 
output signal from 5 and G, respectively. 

The condition (4) in Definition 2 for the system of equations (28) in the set 
Ź(F, m) is equivalent to the condition 


1 10 sg z4(0) J, 
(ra t0)) arm) 


IGZZÓB 


where z,(w) is the first coefficient of Fourier series expansion of the signal z, (t). 
Let us suppose that the first term in (48) attains at w, and the second one — 
at w,. Then we can write 


sup sup ( 


zeZ(F,m) weR 


z,(0) 
al) esi” 


l 20 || 24(04) | 
s ( ( Flo) to) F(jo,) 
l z4(07) 
W> - 4 49 
-||rgaz* WB F(io,) (49) 
Let 
z,(0) a 
oe |FGO)| 


7 — Szczecińskie Roczniki Naukowe t. IX, z. 1 
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for some z and F. Then, there exist w; and m, such that 


z,(04) = 
A F(jo4) 
and 
l 21(04) Si 1 ) || z1(0,) 
+ =mMmy + )4(04) . (30) 
> We ») Fio,| ""MFde)**'*") |IFde,| | 
There also exist w» and m, such that 
z1(02) m 
Ę F (j0)| 
and 
1 | Z1(02) ( 1 | 2,(02) 
+ 0) — =m —— +) (02) r ; (51) 
UE ANA ») F (jo|  *NFóo) © 72)  ||FGo;) 
Substituting (50) and (51) into the inequality (49) we obtain 
1 m]. 1 SA 
su m + )ą(0 —m|| ———+ ),(03) )<o 
sw | nz Z ») (ada A a) 
Or 
i sai l PA 
su —+ )4(0 ) - (sqaz* p, (w: ) )< 0. 
AIG F(jo i0,) Ji : F(j02) A * 
This inequality is equivalent to the inequality 
1 A Ah 1 =" 
su ——— + 71(0;) — su +hloy js 52) 
„wp ((rqoj Ja : ZA (Gio 2) ; E 


Introducing 0;, w, insted of o; , w; and W,, 04 we establish that the inequali- 
ty (52) takes the form (47). Thus we have proved equivalence of conditions (46) 
and (4) in the approximate system (28) for z(t)eŹ(F, m). 

Repeating the same arguments as in the proof of Theorem 1 one can show 
equivalence of conditions (46) and (4) for the system of equations (7) in the class 
of signals Ź(F, m). This consideration can be generalized onto the case when 
Z(F, m) is a set of almost periodic signals. 


EXAMPLES AND CONCLUDING REMARKS 


It follows from Theorems 1 and 2 that the optimal controller and strongly 
optimal one for the system described by the approximate equations 


d y(t) = G(x (0) 
(53) 
x(t) = z(t) — fyl(t-2af (z) 

0 
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are also optimal and strongly optimal, respectively, for the system (7). 
The possibility of replacement G by G significantly simplifies the optimization 
procedure for the considered control system. 

It should be noticed that the determinant in formulas (27) and (46) is not- 
hing else but the Jacobian of operation ©. 

The proofs of Theorems 1 and 2 provide a simple geometric interpretation 
of formulas (27) and (46). 

In the case of minimization of the performance index [,(G) given by (46), for 
every G we take into account the least distance between the curve 1/F(jo) and 
family of curves — G(jo, N). The least distance does not need to be taken for 
the same w in both curves. From controllers Ge C we choose a controller Gy 
for which the least distance between curves 1/F(jw) and —G(jo, N) is a great 
as possible. 

Let us take into account the performance index f,(G) given by (27). 
We look for the least distance between the curve 1/F(jo) and the familly of 
curves — G(jw, N). In this case the distance is measured for the same w in both 
curves 1/F(jo) and —G(jo, N). It may happen that curves 1/F(jw) and 
— G(jo, N) cross for different w. From controllers Ge C we choose a controller 
G, for with the least distance, in the above sense, is as great as possible. 


Example 


Let the plant be described by the transfer function 


l 


As C we assume the set of controllers described by the equations 


y(t) = a(o (0) + BSE A>0, B>0, (55) 


where 


—m*+3m(x+m) for x<-m 
GO" for _m<x<m (56) 
m? + 3m* (x —m) for x>m. 


The function g(x) is continuously differentiable. The describing function of the 
controller GeC takes the form (for x(t) such that |x(t)| < m) 


G(jo, N) = AN? +jABo. (57) 


The class Z(F, m) is defined by the formula 
sup |z(s)(1+sT)| < m 


res>0 
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-Gz63G,6, | 1+ jw 
Fig. 4. Spectrum of operation F" '(s) = 1+ST 
and family of curves — G(ja, N) 


Figure 4 shows the spectrum of operation given by (54) and a family of 
curves — G(jw, N) for various N. Parameters of the optimal controller as well 
as of the strongly optimal one can be evaluated from the formulas 


k E 
I,(G) = inf |+loT+; AN +JABTe 
ŚR 
, | SZM] 
I,(G)= "inf |-+jo, TrzAN +jABTO,|. 
01,02ER 
N 


ADAM ŁOZOWICKI 


O OPTYMALIZACJI NIELINIOWYCH UKŁADÓW ŚLEDZĄCYCH 
W PRZESTRZENIACH H, 


Streszczenie 


W pracy rozważa się nieliniowe układy śledzące. Podano definicję optymalności i mocnej 
optymalności układów śledzących. Zakłada się, że klasa sygnałów wejściowych (zadanych) należy 
do przestrzeni H,. Krótko omówiono metodę funkcji opisującej. Następnie zostały zaprezentowa- 
ne dwa twierdzenia umożliwiające znalezienie optymalnych (w podanym wyżej sensie) nieliniowych 
regulatorów. 


AJIAM JIO3OBMLIKA 


O56 ONTHMAJIMJALIMM HEJIMHEMHbIX CJIENAIIHMX CHCTEM 
B IIPOCTPAHCTBAX H, 


Pe3roMme 


B pa6oTe paccMaTpHBałoTca HeJIMHeHiHbie CJIeNALIHE CHCTEMbl. IIpeHCTaBJIeHbi Oornpene- 
JIEHHA ONTUMAJIBHOCTH MH CHJIBHOŃ ONTUMAJIBHOCTH CJIEĄAINMX CHCTEM. IIpermoJaraeTca, uTo 
KJIACC BXOJNHbIX CHTHAJIOB CHCTEMbl IIDHHAJIJIEXKUHT K IIDOCTpaHCcTBy H,. KpaTko paccMOTpeHbl 
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METOJIbl OTIHCHIBAFOLIEH (byHKLUIMM. KpoMe Toro HpeuIoxeHbl JIBe TeOpEeMbi, IPH IIOMOLINH 
KOTODBIX MOXHO HaliTM OITUMAJIbHbIŃ M CHJIEHO OITHUMAJIbHbBIŃ (B yKA3AHHOM CMBICJie) 
HeJIMHeHHbIH peryJiaTOp. 
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VIBRATION STABILITY OF MACHINE TOOLS 
IN THE PLANNING-DESIGNING PROCESS 
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Summary. The principles of multidimensional machine-tool system vibration 
stability calculations are presented in this paper. The vibration stability analysis is 
composed of two steps: 1 — determination of the absolute limit of vibration stability, 
2 — detection of so called *weak point” in machine tool structure i.e. element that is 
responsible for the loss of vibration stability. The algorithm for determining the limit 
cutting parameters and the complex form of the machine-tool system vibration at the 
stability limit provides a base for developing a computer calculation system that 
assists the machine-tool planning-designing process. 


Key words: machine-tool — vibration stability — cutting process. 


INTRODUCTION 


In the 80s, an experimental method for analyzing vibration stability of 
machine tools was developed and implemented in the research practice at the 
Institute of Mechanical Technology of the Technical University of Szczecin [4]. 
The essence of this method consisted in that the prototype of a machine tool 
was put to an unstable state which appears in the form of self-excited vibration. 
Under this state, vibration was measured for particular solids (considered as 
rigid solids) of a mass-dissipative-spring system (the MDS system) of machine 
tool. The results of measurements were processed using a computer, and then 
the computer animated film was generated. On this basis it was possible to 
evaluate qualitatively the dynamic properties of the machine tool. An additio- 
nal analysis made on the effect of vibration of particular solids on the vibration 
of tool to a workpiece enabled quantitative evaluation and final identification 
of that element of machine tool — chuck — workpiece — tool (MCWT) system 
which was responsible for the loss of vibration stability (so called *weak 
point”). A modification made to the design of this "weak point” contributed 
to the improvement in the vibration stability of the machine. The method has 
been appreciated by machine tools manufacturers. However, in the course of 
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the investigations, a question, posed frequently by the designers, has arisen: 
couldnt such investigations, even though estimative, be made during the 
planning-designing process, since the necessity to build the prototype of machi- 
ne tool for investigation purposes contributes considerably to the increase in 
the project cost. 

To answer such a posed question, it is necessary to posses reliable mat- 
hematical models that describe both the MDS system of machine tool and 
the processes which occur in it — especially the cutting process. Moreover 
these models have to be written in a uniform convention. At the Institute 
of Mechanical Technology of the Technical University of Szczecin a method 
for modelling the MDS system was developed. The principles of this method 
were based on the Rigid Finite Elements Method (the RFE Method) [3] and 
on an original, our own method to model contact joints of slidewise [6]. 
It was, therefore, to: 

— develop principles for modelling the cutting process in multidimensional 
dynamic MCWT System, and 

— develop principles for analyzing vibrational stability of multidimensio- 
nal MCWT System in the planning-designing process. 

Realization of so stated targets is the subject-matter of this paper. 


MODELLING CUTTING PROCESS 
IN MULTIDIMENSIONAL DYNAMICAL MCWT SYSTEM 


The essence of modelling the cutting process in the multidimensional 
MCWT System consists in that the mathematical model of cutting process 
dynamics, which is defined in a coordinate system 0)'x41x2x3 (Fig. 1b) connected 
with the apparent point of contact between a tool and a workpiece, is deter- 
mined in generalized coordinates system, 1.e. those in which the mathematical 
model of the MDS System is defined (Fig. la). Distinguished in Fig. la ( by 
giving the full description) are: the reference system 0„x;,x2,x3, of a RFE No 
r that represents a cutting tool and the reference system 0,Xx15x2,x3, of RFE 
No p that represents a workpiece. 

The transition from one reference system to another is effected by means 
of a transformation matrix. The principles of forming such matrices are 
presented in the work [7]. Figure 2 shows the block diagram of dynamic 
MCWT System which reveals a coupling between the MDS System and the 
cutting process. The generalized displacements u are transformed, by means 
of the transformation matrices N,,, from the generalized coordinates referen- 
ce systems to a reference system that is coupled to the apparent point 
of contact between a tool and workpiece. These displacements (u,„,) produce 
changes in given cutting parameters y, which causes changes to the cutting 
force components P'. The cutting force components are transformed to the 
generalized force Q,„ which produces displacements u in the MDS System. 
The mathematical model of the MDS System is represented by the block 
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r - RFE= cutting tool + gripper 
p-RFE = workpiece 


Fig. 1. Reference systems in which the following mathematical models are defined: a) machine tool 
MDS System model; b) cutting process model 


(MDS) 
MDS-MCWT System 


Wyps 


Matrixthattransforms 
forces from tool system 
to generalized 
displacements system 


Matrix that transforms 
generalized displacements 
to tool system 


i 
Ni 


Fig. 2. Block diagram of dynamic MCWT System with cutting process feedback 
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Wups — a matrix of dynamic characteristics of the MDS System which are 
determined being based on the RFE Method [3]: 


Wyps = (M:s*+H's+K)"*, (1) 
where: 


M = mass matrix, 

H = dumping matrix, 

K = stiffness matrix, 

s = Fouriers transformation operator: s = jw, where w = circular frequen- 
cy, j=1maginary unit. 

The cutting process model is represented by the block Wys (matrix of dyna- 
mic characteristics of cutting process [5]). Such a way to model the cutting 
process enables many *coupling elements” in form of cutting process is to be 
considered in the vibration stability analysis. Generally, many cutting processes 
can be executed within the same dynamical MCWT System, which results in 
certain implications for the presented in block diagram (Fig. 2): every new 
cutting process creates a new feedback loop with a new transformation matrix 
N.,, and, moreover, introduces cutting parameters p and matrices of dynamic 
characteristics of cutting process Wpę and finally, causes changes to the force P' 
components. This affords an opportunity for a multi-variant analysis of 
vibration stability of multi-tool and multi spindle systems and, what is more, 
there is a possibility to carry out such an analysis for the sake of one chosen 
cutting process, all the processes simultaneously or certain combination of 
them. 

There remains, however, the key problem of inner structure of a block that 
represents definite cutting process. This structure depends on the process exe- 
cuted on the machine tool. Below are presented examples of forming this 
structure for two selected cases of cutting processes: single-tool turning and 
multi-tool turning. The rules assumed here may be extrapolated to other cases 
of machining. 


SINGLE TOOL TURNING 


The scheme of constrained cutting zone at single tool machining is shown 
in Fig. 3. 

The relationship between changes in the cutting force components and 
displacements of the tool relative to the workpiece is described by the following 
formula: 


P'(s) = Wz(5): Ly (5) — ur(5)], (2) 


where P' (without the Fouriers transformation operator symbol — which 
is omitted in subsequent designations) determined as: 


P' = coliP4, Pa, P3, 0, 0, 0; (3) 


Cutting process 
(inner undulalion) 


Cutting process 
(outer undulalion) 


Fig. 3. Block diagram of cutting process at single-tool turning 


is the vector of cutting force components, and: 
l ' ! ! ' ' I 
Uw r= col (Uefa Uwr2> twyas Uwr4; Uwr5; Uwr6j (4) 


is the vector of displacements of tool relative to workpiece, whereas: 


WyR11 WyR12 WyR13 0 0 0 
WpR21 WóyR22 WóyR23 0 0 0 
WyR3 1 WyR3 2 WyR3 3 0 0 0 
R = 5 
Hn=|lo 0 6 600 (5) 
0 =0 %05 6% 06 
0.0.0 000 


is the matrix of dynamic characteristics of cutting process. 

The vector y =colity;, V>, Va, V4> Vs, VG) designates the given values for 
reciprocal position of tool and workpiece to provide realization of cutting 
process. 

It is to be stressed that the relative displacements u,, influences the cutting 
force components in an indirect way, namely, the displacements produce 
changes only in the cutting parameters (thickness of cut a, width of cut b, 
and cutting speed v) and these, on the other hand, cause changes to the 
cutting force components. The cutting parameters depend on displacements 
actually occurring at a chosen moment of time t (so-called inner undulation 
u„,(t)) and on what is the irregularity of cut surface that was caused 
by displacements occurred at the time (t—T), i.e. one revolution earlier 
(so-called outer undulation u,„,(t—T). Simultaneous occurrence of both 
forms of changes in the cutting parameters is termed a vibration-ridged 
work —surface machining (vibration —ridge regeneration). The cutting force 
components are therefore a sum of changes resulting from both the undula- 
tions: 


P' = Pp+ Ppr (6) 
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and, at the same time: 
P, = Wp'(y —u,,), (7) 
pr = Wpr'lw,'€ *', (8) 


where: 

W, Wpr = matrices of dynamical characteristics of cutting process that re- 
late to inner and outer undulations of the thickness of cut respectively; the 
structure of these matrices is identical with that of matrix Wpę, but, for the sake 
of different forms of changes in the cutting parameters, they may differ from 
each other. 

Presented above relations between the cutting force components and the 
relative displacements of tool and workpiece can be represented by means of 
a block diagram, as shown in Fig. 3. 

From comparison between Fig. 3 and Fig. 2 it can be concluded that: 


Wpę = —Wp+ Wpr'e 77. (9) 


One of more important aspects in forming the matrices Wp and Wpr is to 
select the form of the mathematical relationships which describe particular 
terms of the matrices Wp and Wpr — the so-called dynamical characteristics of 
cutting process (DCCP). In the world specialistic literature there is no explicit 
opinion on this subject and the physical interpretation of dynamic phenomena 
occurring in the cutting process as well as the resulting final relationships are 
the subject of scientific disputes [2]. In the work [7] an analysis is made for 
several extremely different mathematical models of cutting process which were 
constructed based on various DCCP. It seems that the directives, suggested in 
that work, for selecting the forms of relationships that describe particular terms 
of the matrices Wp and Wpr may be useful in analyzing dynamical properties of 
machine tools in the planning — designing process. The directives recommen- 
ded to assume the so called [7] proportional model of cutting process for 
computations of the stability. 


MULTI-TOOL TURNING 


The multi-tool turning process is characterized by that, at every moment, 
cutting is performed by several cutting edges simultaneously (Fig. 4). This 
complicates the task to model the cutting process dynamics for this case of 
machining. 

By developing the modelling concept assumed for single cutting edge so as 
to cover the multi-tool turning instance the relations between the relative disp- 
lacements vector u,, and the force vector P' can be represented by means of 
a block diagram shown in Fig. 5. Both the vector u,, and the vector P' are 
defined in selected reference system (' x; x+x3; (in Fig. 4 this system has been 
coupled with the apparent point of contact between the tool No 1 and the 
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Gripper (RFE no r ) 


Fig. 4. Multi-tool turning 


workpiece). The displacements u,,. are transformed to the system coupled with 
the apparent points of contact O; x;4 x; x;3 of particular cutting edges with the 
workpiece. This transformation is effected by means of the transforming mat- 
rices G„,;. The cutting forces exerted on particular cutting edges are transfor- 
med by means of matrices G;,,, to the reference system 0'x; x+x3 and summed 
there. In this model the effect of vibration-ridged work surface machining is 
taken into account and, moreover, the time constant T is identical. Based 
on the presented block diagram it is possible to derive a relationship which 
represents the mathematical model of cutting process dynamics for this type of 
machining: 


m m 
/ , —ST T ==$'JF T 
Wy = —Węt Wrce "=—) Gui WiGuite *"* 2, Gui WeriGwri 


i=l1 i=1 


= > Gui: (—WpitWprs'e *')'Gi, (10) 
i=1 
where: 

Wp; = the DCCP matrix related to the inner undulation and expressed in 
the reference system coupled with the apparent point of contact of i-th cutting 
edge with the workpiece (6 x 6); 

Wqr: = the DCCP matrix related to the outer undulation and expressed in 
the reference system coupled with the apparent point of contact of i-th edge 
with the workpiece (6 x 6); 

m = number of cutting edges engaged simultaneously. 

The rules of constructing the above matrices are the same as for single tool 
turning. The principles of constructing the matrix G,,, are presented in the 
work [7]. 
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Fig. 5. Block diagram of cutting process model for multi-tool turning 


ANALYSIS OF MCWT SYSTEM VIBRATION STABILITY 


In the planning-designing process made on a machine tool, the vibration 
stability of a closed multidimensional dynamic MCWT System (Fig. 2) is ana- 
lyzed based on the Nyquist stability criterion which is expressed by the rela- 
tionship: 


det[T' + Wyps' Wea] = 0 (11) 
Or 
det[I+ Wyps: WpęR] = 0, (12) 
where: 
I,I =unit matrices, and 
Wyups = Ni: Wiups” N wr» (13) 
WR = N,,' Wpę' Nur. (14) 


The equations (11) and (12) are identical to one another but differ in 
the order of determinant: the complex determinant (11) is of order 6 whereas 
the complex determinant (12) is of order 6n, where n is the number of RFEs 
in the MDS System model. The dimension of the matrix N,, is 6nx6. 
This fact is essential for the vibration stability analysis which is composed 
of two steps. 
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Step 1 


This step consists in determining the absolute limit of vibration stability, 
1.e. determining a surface of limit widths bg, of layers cut by particular 
edges: 


J (bo1, boa» Bos; ---, bom) 5 0 (15) 


such that enables the equations (11) and (12) to be satisfied at any combination 
of the widths of cuts by, from this surface. Moreover, it is required that the 
criterion for the minimum sum of widths of layers cut by particular edges to be 
also satisfied: 


y = 3 bo: = min. (16) 
i=l1 
Particular widths bg, are related to the equation (11) in the following way: 
after substituting the equation (10) to the equation (11) we obtain: 


det [T' + Wyps 2, Gui: (— Wpi+ Wprs:e *'): G„„| = 0. (17) 
i=l1 
Since the characteristics Wp; and Wpr, are proportional to the width of cuts bg;, 
we may write down that: 


" 24 
W; = bo: " WPi» 
Wprq: = Doi' WpTi- (18) 


By substituting the (18) to the (17) we obtain: 
det [r + Wyns 2 GT; ; ( = Wp; aj WPT: | 8,4 j Gui i boi | =0 (19) 
i=l1 


The determination of the surface (15) is therefore reduced to solving a clas- 
sical optimization task in which, at preassumed combination of values 
boi (i=1,2,...,1—1,1+1,..., m), such that it provides vibration stability of 
MCWT System, the remaining value of by, is selected so that the system is 
brought to the limit of vibration stability. The essence of this search consists in 
minimization of the function (16) at the equality constraints (19) and with the 
inequality constraints being retained 


dnia s, T < diz: 
(20) 
bo; 2 0, 


where: 

C Tmax» Tmin> = spindle revolution period interval which should satisfy the 
condition f'(Tnax — Tin) 5 1, where f is designated that vibration frequency 
which is being currently considered. 
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The entire calculation cycle is repeated for successive values of frequency 
J that correspond to the table of discretely determined matrices of frequency 
characteristics Wyps. For the calculation that frequency and matrix are selected 
at which the sought value of bg, is the least. Next, the same cycle is repeated for 
other combination of by,. Each combination allows individual point to be 
determined on the surface (15). One of the fundamental problems which arise in 
realization is the calculation time especially the time to calculate the complex 
determinant (19) which regardless of the optimization method assumed, must 
be determined in every iteration of the optimization procedure. It would there- 
fore be advisable that this determinant is of possibly low order. This is why it is 
convenient, at this stage, of vibration stability limit determination to make use 
of the formula for the condition of vibration stability written in the form (11) 
and in its detailed expansion — the form (19). 

Realization of the step 1 allows the designer to answer the question: whet- 
her the newly designed machine tool is resistant to self-excited vibration (the 
loss of vibration stability) within the foreesen ranges of cutting parameters. 
If the answer is positive, the vibration stability analysis is terminated at this 
step. If the verdict is negative, then it is needed to modify the structure of newly 
designed machine tool in order to improve its vibration stability. This modifi- 
cation can be introduced based on extended information obtained in the step 2. 


Step 2 


This step consists in detection of the so called *weak point” ie. an element 
that is responsible for the loss of vibration stability in the MDS—MCWT 
System. The modification made to the structure of this element should cont- 
ribute to the improvement of dynamical properties of the entire MDS System. 
In order to detect the weak point you should select any, known from the 
predicting calculation step, combination of values bg, from the surface (15), so 
the one of which satisfies the determinant (19) nulling condition, calculate the 
matrix Wp; and Wpr; at the vibration stability limit (the formula (18)) and, next, 
based on the formula (10), calculate the matrix Wpę at the vibration stability 
limit for the period T resulting from the loss-of-stability condition. 

Having determined the matrix Wpą it is possible, on the base of the formula 
(14) to determine the matrix Wpę. The matrix satisfies the determinant (12) 
nulling condition. By multiplying this determinant by the set [Wyps] we 
obtain: 


det[ Wyps]' det[I + Wyps* Wea] = 0, (21) 
which means that: 
det[Wyups + Wre] = 0, (22) 
where (cp. formula 1) 


Wyjs = M:'s?+H's+K. (23) 
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The square brackets in the equation (22) embrace a matrix of coefficients of 
system of homogeneous complex equations. This system of equations describes 
vibration of the system u at the vibration stability limit: 


(M:s+H's+K+ Wyp): u=0 (24) 


Therefore by solving this system of equation we obtain a complex form of 
vibration u at the vibration stability limit and moreover: 


U, a Ay; i Mb: (25) 


where: 

Au;, ©;, = amplitude and phase of vibration respectively in the direction 
i-th generalized coordinate. 

Vibration of the system in the time domain are described by the following 
system of equations: 


u,(t) = Ai" SIn(Ot + Qu), (26) 


where w is circular frequency of self excited vibration. 

Based on the equations (25) and (26) it is possible to: 

— generate a computer animated film that illustrates the image of machine 
tool vibration in terms of machining at the vibration stability limit, 

— indicate those elements of the MCWT System which have the strongest 
effect on the amplitudes of relate vibration of tool and workpiece, that is the 
„weak points” in the MCWT System. 

The above tasks can be realized by means of a computer program 
„MODAN” [1] which is available in the Institute of Mechanical Technology 
at the Technical University of Szczecin. 


CONCLUSIONS 


The presented method for computing the vibration stability of a multi- 
dimensional MCWT System enables, apart from the vibration stability estima- 
tion issued in the form of limit widths of cuts, the determination of complex 
forms of vibration at the vibration stability limit. Hence this enables the „weak 
point” to be indicated in the MCWT System. This information is of great 
practical significance, since it facilitates the decision to be taken, by the desig- 
ner on possible modification made to the developed structure just in course of 
the planning — designing process. 

The algorithm proposed in this work for determining the limit cutting 
parameters and the complex forms of the MCWT — System vibration at 
stability limit makes a basis for developing a computer calculation system that 
assists the machine tool planning-designing process. 
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JAN TOMKÓW 


TEORETYCZNE PODSTAWY PROGNOZOWANIA 
WIBROSTABILNOŚCI OBRABIAREK 
W PROCESIE PROJEKTOWO-KONSTRUKCYJNYM 


Streszczenie 


W pracy przedstawiono zasady obliczeń wibrostabilności wielowymiarowego układu obra- 
biarka-uchwyt-przedmiot obrabiany-narzędzie. Analiza wibrostabilności składa się z dwu etapów: 
1) wyznaczanie absolutnej granicy wibrostabilności obrabiarki; 2) wykrywanie tzw. „słabego 
ogniwa” w strukturze obrabiarki, tzn. elementu odpowiedzialnego za utratę wibrostabilności. 
Algorytm wyznaczania granicznych parametrów skrawania oraz zespolonych postaci drgań na 
granicy wibrostabilności tworzy podstawy do opracowania systemu programów komputerowych 
do wspomagania procesu projektowo-konstrukcyjnego. 


AH TOMKYB 


TEOPETMUECKME OCHOBBI IIPOTHO3MPOBAHNHA 
BHBPOYCTOKUMBOCTM CTAHKOB 
B IIPOEKTHO-KOHCTPYKIIMOHHOM IIPOLIECCE 


Pe3oMe 


B paóoTe rnper1cTaBJIeHbl HPHHIIMIIbI BBIUHCJIEHMA BHOPpOyCTOŃUMBOCTH MHOTOMEpHOŃ CHC- 
TEMbl CTAHOK-IDPEJHMET-HHCTPyMEHT-JIepXKaBKA. AHaJIH3 BHODOyCTOŃUMBOCTH COCTOHT M3 ĄBYX 
3TANOB: 1) onpenejleHne aOCOJIIOTHOTO HNpenejla BHOpoycTOŃUMBOCTH; 2) BbIABJIEHKE, TAK HA3bIBA- 
eMo, „,cJIaOofńi TOUKku” B CTpykType CTAHKa, T.€. 3JIEMEHTA, H3-3A KOTODOTO BO3HHKAET IIOTEpA 
BHOpOycTOŃUMBOCTH. AJITODHTM OIIpEJEJIEHHA IIpETEJIbHbIX IIApAMETPOB pe3AHHA H KOMIIJIEKCHbIX 
BUJIOB BHUOpAaIHM CO3JIA€T OCHOBbI JIJIA pa3paQOTKH CHCTEMbl KOMIIbIOTEDHBIX IIDOTPAMM JIA 
NONIEpKKM IHPOeKTHO-KOHCTDYKIIAOHHOTO IDpOLECCA. 
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